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On  Direct  Methoda  for  Solving  Syrnmetrlc  Systems 
of  Linear  Eguatlona 
by 

James  Raymond  Runoh 

Abstract 

There  has  been  no  stable  direct  method  for  solving  symmetric 
systems  of  linear  equations  wtilch  takes  advantage  of  the  symmetry.  If 
the  system  Is  also  positive  definite,  then  fast,  stable  direct  methods 
(e.g.,  Cholesky  and  synanetrlc  Gaussian  elimination)  exist  which  preserve 
the  symmetry.  These  methods  are  unstable  for  symmetric  indefinite  sys- 
te  IS.  Such  systems  often  occur  in  the  calculation  of  eigenvectors. 
Gaussian  eiiminatlun  with  partial  or  complete  pivoting  is  currently 
recommended  for  solving  symmetric  indefinite  systema,  and  here  symmetry 
is  lost . 

We  present  a  generalization  of  symmetric  Gaussian  elimination, 

called  the  diagonal  pivoting  method,  in  which  pivots  of  order  two  as 

well  as  one  are  allowed  in  the  decomposition.  We  show- that  the  diagonal 

pivoting  method  for  symmetric  indefinite  matrices  takes  advantage  of 

symmetry  to  that  only  ^  n^  laultipllcations ,  at  most  -j  n^  additions, 

1  2 

and  n  storage  locations  are  required  to  solve  A  x  -  b  ,  where  A 

is  a  non-singular  synmetric  matrix  of  order  n  .  Furthermore,  w«  -show 

^  > 

Chut  Che  method  is  nearly  as  stable  as  Gaussian  elimination  with  complete 
pivoting,  while  requiring  only  half  the  number  of  operations  and  half 
the  storage. 


We  Include  a  listing  of  an  Algol  procedure  for  Che  diagonal 
pivoting  nethodt  which  Is  applicable  both  Co  sytninetrlc  definite  and 
Indefinite  systems. 

We  discuss  the  problem  of  symmetric  band  matrices  and  present 
an  algorithm  only  for  Che  crldlagonal  case.  Further,  we  discuss  Che 
problem  of  equilibrating  symmetric  uatrlces  while  preserving  symmetry 
and  we  present  a  simple  algorlttus  (and  Algol  procedure)  for  accomplishing 


this . 
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Chapter  1  ;  Introduction 

1 . 1  Syinm'»trlc  Syatema  of  Linear  Lquatlona 

Let  ua  conalder  direct  methods  for  solving  the  system  of  linear 
algebraic  equations,  A  x  ~  b,  where  A  is  symmetric. 

If  A  is  aiao  positive  definite,  then  Cholesky's  method  (§2.7) 
and  L  D  L^  method  (§2.6)  are  fast,  stable,  and  preserve  symmetry. 

T.f  A  Is  symmetric  but  indefinite  (neither  positive  definite  nor 
negative  definite),  Cholesky's  method  and  the  L  D  L**  method  are  unstable 
end  can  produce  very  Inaccurate  results  (§2.8). 

At  the  present  time,  if  A  is  symmetric  indefinite,  Gaussian 
elimination  with  partial  or  complete  pivoting  Is  recommended  for  solving 
the  system  (Fox,  p.  80,  185),  and  thus  the  symmetry  of  A  is  ot  no 
advantage . 

Is  there  an  algorithm  for  the  symmetric  Indefinite  case  wh ' ch  is 
stable,  is  faster  than  Gaussian  elimioation,  and  can  take  advantage 
of  the  symmetry? 

1 . 2  Our  Contribution 

We  diacuas  the  problem  is  Chapter  2  and  review  previous  efforts  in 
Chapter  3.  In  Chapters  4-6  we  present  a  method,  called  diagonal  pivoting, 
which  fulfills  the  above  requirements  when  reetrlcted  to  equilibrated 
natrlces .  In  Chapter  7  we  present  a  method  for  equilibrating  symmetric 
matrices  In  a  very  simple  manner.  A  variation  of  the  dia(,onal  pivoting 
method  is  presented  in  Chapter  tij  this  method  is  applicable  to  unequlll- 


2 


braced  laacrlces  and  it  fulfllla  the  above-mentioned  requirements.  In 
Chapter  9  we  show  that  the  diagonaL  pivoting  raechod  is  almost  as  fast 
as  Cholesky  or  L  D  .  In  Chapter  10  we  perform  a  backwards  error 
analysis.  In  Chapters  11-12  we  show  that  the  method  is  essentially  as 
stable  as  Gaussian  elimination  with  complete  pivoting  (in  the  sense  of 
Wilkinson's  analysis  for  Gaussian  elimination  with  completely  pivoting, 
Wilkinson  (1961)).  In  Chapter  13  we  show  that  iterative  Improvement  is 
as  applicable  here  as  it  is  for  Gaussian  elimination .  In  Chapter  14 
we  discuss  the  problem  of  symmetric  band  matrices. 

All  the  results  proved  are  applicable  to  complex  systems  where  A 
is  Hermitian. 

1.3  Summary 


’.et 


A  be  an 


n  ^  n  matrix  with 


max  |a  I 
ij 


1  . 


We  want  to  solve  A  x  ■»  b  . 

k.  i 

Let  ~  C,  denote  C,  N  ^  C,  N  ,  where  the  C.  ,  0  ^  i  £ 

1-0  ^ 

ere  constants  independent  of  n  . 

Let  G.E.  denote  Gaussian  nliminatlcn. 


The  situation  is  sunznarlzed  in  the  following  table: 


- - 

Method 

r- 

Reatrlc— 

tio'.ia 

on  A 

Number  of 

Mult j  olications 

Number  of 
Additions 

Storage 

Bound  on 
Element 
Growth 
(Stabi llty) 

G.E.  with 

complete 

pivoting 

det  A  0 

~  n* 

f  (n) 

G.E.  with 

partial 

pivoting 

det  A  0 

Cholesky 

Method 

symmetric 

positive 

definite 

“i"’ 

~  1  _2 

2  " 

1 

Diagonal 

Pivoting 

symmetric, 
det  A  0 

'h‘ 

~  1  n2 

2  " 

/n  f  (n)  X 
c (a)  h(n,a) 

fn  —1 

Here  0  <  a  <  1  ,  f  (n)  “  >  II  k  ^  * 

lk.2  J 

dependent  on  the  pivotal  strategy,  and  c<a)  and 
§§11.6-7.  In  Chapter  12  we  show  that  c(a)h(n,a) 


h(n,a)  Is  a  function 

h(n,a)  are  defined  in 
<  3.07(n-l)°‘^'^^  <  2 


for  a  -  (1  +  /r7)/8  -  a  . 

o 

1 . 4  Origin  of  Symmetric  Indefinite  Systena 

The  problem  of  indefinite  systems  cf  linear  equations  is  sometimes 
dismissed  as  academic  by  the  claim  chat  physical  pr.^blems  always  generate 
positive  definite  systems  of  linear  equations.  However,  the  numerical 
solution  of  natural  problems  often  gives  rise  to  situations  which  do  not 
have  a  physical  origin.  We  give  two  related  examples. 


1 


t 


In  the  Rayleigh  Quotient  Iteration  for  finding  eigenvalues  of  a 
positive  definite  matrix  A  CWllk-inson  (1965),  p.  172,  and  particularly 


1 

3 


4 


p.  629)  we  need  to  aolve  the  systems  (  A  r^I  )  »  where 

★  * 

''1  ”  *1  ^  *1^  **1  *1  '  ^  cannot  be  definite  because 

lies  between  the  extreme  eigenvalues. 

In  the  inverse  Iteration  method  for  finding  the  eigenvector  corres¬ 
ponding  to  an  approximation  X  to  an  intermediate  eigenvalue  of  a 
positive  definite  matrix  A  ,  we  need  to  solve  (  A  -  X  I  )  »  u^  , 

u^_^^  “  '’t+l'^  ^'^i+1^  ■  A  -  X  I  can  be  indefinite  even  when  A  is 


positive  definite  (Wilkinson  (1965),  pp.  618-635). 
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Cliapter  2  :  Presentation  of  the  Problem 

2 . 1  Introduction 

The  speed  and  storage  capacities  of  current  digital  computers 
allow  us  to  solve  large  systems  of  linear  equations  by  direct  methods. 

Here  we  shall  consider  direct  methods  for  solving  a  system  of  linear 
equations,  A  x  *  b  ,  whe'e  A  Is  symmetric  and  det  A  0  . 

2.2  General  Problem  :  Exact  Arithmetic 

First  let  us  consider  the  solution,  in  exact  arithmetic,  of 
A  X  »  b  ,  where  A  Is  general  and  det  A  0  .  We  know  that  Gaussian 
elimination  will  give  ttie  aolutlon  provided  that  whenever  a  zero  appears 
In  Che  leading  diagonal  position  we  Interchange  Chat  row  with  a  lower 
row  with  non-zero  leading  element  (such  a  row  will  exist  since 
det  A  j*  0),  e.g.  if  A^^j^  «  0  and  if  J  Is  the  least  integer  for  which 
Ajj^  ^  0  ,  then  we  Interchange  rows  1  and  j  .  Or,  equivalently,  there 
exists  a  permutation  matrix  P  sucli  that  Gaussian  elimination  without 
interchanges  applied  to  PA  will  give  us  the  solution. 

Since  we  could  also  do  the  same  with  columns  Instead  of  with  rows, 
there  exists  a  permutation  matrix  Q  such  that  Gaussian  elimination 
without  Interchanges  applied  to  A  Q  ’’ill  give  us  the  solution. 

In  matrix  notation  the  Gaussian  elimination  algorithm  factors  A 
Into  A  -  L  U  ,  where  L  Is  unit  lower  triangular,  U  la  upper  triangular 
and  L  and  U  are  unique  when  they  exist  (Wilkinson  (1965),  p.  204). 

Thus  ■--xact  arltlunetic  there  exist  permutation  matrices  P  and  Q 


sucli  that  P  A  »  Lj^  and  A  Q  «  D2  1  provided  only  that 
det  A  1*  0  . 


2 . 3  General  Problem  :  Failure  Previous  Attempts  in  Finite 

Precision  Arithmetic 

In  finite  precision  arithmetic  (Wilkinson,  1963)  the  above  algo- 
ritlrni  can  fall  if  we  interchange  only  under  the  condition  that  the  ele¬ 
ment  in  the  leading  diagonal  position  be  zero. 


Let  A 


Then  A  =  L  U  = 


€  1 

’l  0 
1/C  1 


and  b  “ 

C  1 

Lo  n  -  1/eJ 


1/e 

0 


.  However,  if  e  and  n  are 


small  enougii,  then  in  finite  precision  arithmetic  the  operation  n  -  i/c 
yields  -1/e  . 

Let  U  and  x  be  the  matrix  and  vector  of  the  values  of  U 
c  c 

and  X,  respectively,  computed  in  finite  precision. 


Then  U 


e 

1 

-1 

1/c 

u  » 

c 

0 

-l/e_ 

and  a  ”  L  b  - 

So  X  “  U  ^  a  “ 

‘  -n  " 

-n/c 

e  (l-ne) 

▼ 

1 

l/e 

U(i-nc)-l 

- 

,  but  X  »  U  ^  a 
c  c 


0 

l/£ 


If  n  =  c  c,  then  x  v 


ll/e 


,  and  the  error  in  the  first  component 


of  the  computed  solution  is  Ic]  . 
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2.4  Condition  of  Matrix 

Wilkinson  (1965,  pp .  189-191)  shows  that  the  relative  error  in  the 
solution  of  a  system  of  linear  equations  is  bounded  in  proportion  to 
the  condition  number  of  A  ,  k(A)  ■  IaI  IaI  ^  ^  1  ,  l.e..  If  e  »  x  - 
then  lei  /  Ixl  p(A)  <(A)  .  We  should  not  expect  a  small  error  If 
|>:(A)  la  large. 

-1  1  h 

Here  A  has  a  very  satisfactory  condition  number.  A  ••  -r 


*  “1 

Using  the  one-norm,  IaI  “  max  I  KJ.  we  have  k(A)  ”  “ 

j  1-1 

[1  |_maxCLe|.  |nl)l^  i  , 
fi  -  ne| 


The  computer  replaces  en  -  1  by 


-n  1 

.  Then  (.A  b  « 

-n/e~ 

_  1  -e 

c 

-  1/^. 

and  the  computed  inverse 

.  Thus  the  trouble  lies 


is 

in 


the  Gaussian  elimination  algorithm,  not  in  the  matrix  A  . 


2.5  General  Case  :  Stable  Direct  Methods 

(a)  Direct  Inversion 

Direct  inversion  (the  formation  of  A  ^  )  of  a  system  requires  ~ 
each  of  multiplications  and  additions  (Fox,  pp.  177-179).  Gaussian  eli¬ 
mination,  however,  requires  only  ~  n^  each  of  multiplications  and 
additions.  Thus  we  would  prefer  to  use  Gaussian  elimination  if  we  could 
obtain  a  satisfactory  solution. 
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(b)  Stability 

Let  us  attempt  tu  solve  A  x  b  ,  where  A  is  nxn  and 
det  A  ft  0  .  If  Is  Che  solution  we  obtain  from  the  computer,  we 

may  consider  to  be  the  exact  solution  of  the  system  (A  +  E)  y  “  b  . 

Ue  might  say  Chat  the  algoclchm  we  use  Is  stable  if  the  elements  of  E 
ace  small  In  comparison  to  Che  corresponding  elements  of  A  .  Actually 
the  term  Is  more  often  used  when  IeI  /  lAl  is  small.  (Here  I  I 
Is  any  norm.) 

(c )  Stability  for  Gaussian  Elimination 

Uillclnson  (1960)  showed  Chat  for  Gaussian  elimination  we  have; 

Ie  I  _<  2.01  max(l-l,J)2  ^  max  A^”  ,  where  t  Is  the  number 

'  '  r  ,s 

(n”k+l ) 

of  binary  digits  In  the  machine,  k  “  mln(l,j),  and  A^  is  the 

reduced  matrix  of  order  n-k+1  In  the  elimination  process. 

The  important  lesson  from  the  above  Is  that  we  must  be  interested 
in  keeping  the  elements  in  the  reduced  matrices  small.  There  are  two 
well-known  strategies  for  choosing  permutation  matrices  P  and  Q 
such  that  Gaussian  elimination  without  Interchanges  applied  to  FAQ 
will  provide  sufficiently  small  element  growtii  In  the  reduced  matrices. 

(d)  Complete  Pivoting 

The  first  strategy,  called  complete  pivoting,  requires  that  we 
bring  the  largest  element  in  Che  reduced  matrix  Into  the  leading  diagonal 
position.  This  strategy  is  called  complete  since  we  search  the  entire 
reduced  matrix.  Wilkinson  showed  that  this  complete  strategy  gives 
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max 

i.J 

In  words, 
large;  so 
bound  is 


k.+l)j  ^  K'^f  (Vt)inax 
^  '  i.J 


'IJI 


,  where  £(k.)' 


k 
n  r 
r-2 


1_ 

r-1 


the  elements  in  the  reduced  matrices  can  never  become  too 


tilts  strategy  is  never  bad.  It  is  conjectured  that  the  true 


max 

I.J 


(n-k+l^l 

^J 


<  k  where  A  is  real  (l>12.4). 


Equivalently,  the  above  says  that  there  exist  permutation  matrices 
P  and  Q  such  that  Gaussian  elimination  without  Interchanges  applied 

to  P  A  Q  gives  max  |e^j|  £  2.01  n^^^  £ (n)  2^*'  max  |A^J  . 


i.J 


IJ 


I.J 

(e)  Partial  Pivoting 

The  second  strategy,  called  partial  pivoting,  requires  that  we 
bring  the  largest  element  in  Che  £lrst  column  o£  the  reduced  matrix 
into  the  leading  diagonal  position.  This  strategy  Is  called  partial 
since  we  search  only  a  part  of  the  reduced  matrix.  This  is  equivalent 
Co  Che  application  of  Gaussian  elimination  without  interchages  to  PA, 
where  P  is  a  permutation  matrix.  Here 


max 

i.J 


!A^:‘‘'"^''i<  2‘^ 


IJ 


This  bound  Is  sharp  since  A 


10  0 

-110 
-1  -1  1 


-1  -1  -1 


0  1 
0  1 
0  1 


-1  1 


,  where  A  is 


nxn  ,  has  A^"^  -  2“  .  Thus  max  |e_  |  <  2.01  n  2"'  2"  max  |A,J  , 

""  I.J  ■  1,J 


and  this  is  very  weak  when  n  >  t  . 

Correspondingly,  we  could  use  a  partial  pivoting  strategy  in  which 
we  bring  the  largest  element  In  the  first  row  of  the  reduced  matrix  into 


I 

] 

Che  leading  diagonal  position.  Thus  there  exists  a  permutation  matrix 

such  that  Gaussian  ellrolnaLlon  without  interchanges  applied  to  A  Q 

lias  an  error  matrix  E  with  max  .|  f.  2.01  n  2  ^  2*'  max  |A,  .|  . 

l.j  ^  i.J  ^  I 

(£)  The  Error  Matrix 

Uc  see  Chat  the  error  matrix  E  is  dependent  on  Che  decomposition 
L,U  .  Che  matrix  A  ,  the  right  hand  aide  b  ,  and  the  permutations  P 
and  Q  by  which  we  can  pre-  and  poet-mulclply  A  ,  l.e.  E  - 
E(L,  U,  A,  b,  1‘,  y)  ,  wiiere  L  0  =  P  A  Q  .  (For  the  partial  pivoting 
strategy  on  the  first  column  of  the  reduced  matrix,  we  cake  Q  ■  I  in 
the  above . ) 

2 . 6  Symmetric  Case  :  Direct  Methods 

If  A  is  symmetric,  then  we  can  only  apply  congruences  to  A  if 
we  want  to  preserve  symmetry.  In  particular,  whenever  we  interchange 
two  rows,  we  must  also  Interchange  the  corresponding  columns.  Thus 
only  a  diagonal  element  can  be  brought  Into  the  leading  diagonal  position. 

The  :.ymmetric  form  of  the  Gaussian  elimination  decomposition  L  U  gives 
the  decomposition  L  U  ,  where  L  is  unit  lower  triangular,  D  Is 

diagonal,  and  is  the  transpose  of  L  .  Since  we  may  only  perform 

congruences  on  A  ,  the  error  matrix  Is  E  «  E(L,  D,  A,  b,  N,  N^)  where 
N  is  a  permutation  matrix  auch  that  L  D  -  N  A  . 

2 . 7  Symmetric  Positive  Def Inlte  Case 
(a)  Cholesky's  Method 

The  Cholesky  decomposition  (Wilkinson  (1965),  pp.  229-232)  is  the 


roost  well-known  decomposition  for  a  positive  definite  matrix  A 
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(l.e.  Ax  >  0  for  x  0  ) .  Here  A  «  L  ,  where  L  l8  lower 
crlangulsr.  No  Interchanges  are  re^iulred  for  stability. 

(b)  L  D  Decomposition 

If  A  is  positive  definite,  then  its  L  D  decomposition  (the 
symmetric  form  of  Gaussian  elimination)  is  stable  in  the  absence  of 
underflow  and  overflow.  The  elements  of  L  can  be  arbitrarily  large, 
but  if  they  do  not  overflow  then  in  fact  the  error  matrix  E  is  as 
small  as  the  error  matrix  for  the  Cholesky  decomposition  of  A  . 

(Note:  L  -  L  ) . 

(c)  Method  of  Congruent  Transformations 

This  method  (Westlake,  p.  21;  De  Meersman  and  Schctsmans,  (196A)  )  uses 
decomposition  (b)  with  interchanges.  Here  the  largest  diagonal  element  is 
brought  into  the  leading  diagonal  position  at  each  step.  If  A  is  posi¬ 
tive  definite,  then  the  elements  of  L  are  bounded  by  1  and  the 
method  is  stable. 

(d)  Summary 

If  A  Is  positive  definite,  then  the  above  three  methods  are  stable. 

1  2 

If  A  is  nxn  ,  then  each  method  requires  ~  *2  storage  positions, 

~  n’  multiplications,  and  ~  ^  n*  additions  to  solve  A  x  ■  b  . 

2.8  Symmetric  Case  ;  Failure  of  Cholesky  and  L  D  Methods  in 

Exact  Arithmetic 

If  A  is  symmetric  but  indefinite,  then  the  L  D  decomposition, 
the  method  of  congruent  transformations,  and  the  Cholesky  decomposition 
fail  in  exact  arithmetic  for  a  matrix  as  simple  as 

fo  1 

0. 


A 


1 


,  that  is,  there  exists  no  permutation  matrix  N  such  :hat 


NAN*'  has  an  L  D  L*  or  an  L  L*  decomposlclon. 

Note  that  L  L*  Is  always  positive  aemi-def inlte .  Thus  Cholesky 
decomposition  will  fall  In  exact  arithmetic  whenever  A  Is  Indefinite 
and  det  A  ^  0  ,  l.e.  there  exists  no  permutation  N  such  that  NAN*" 
L  L*  where  L  Is  lower  triangular. 

The  L  1)  L*  decomposition  and  the  method  of  congruent  transfor¬ 
mations  will  full  in  exact  arithmetic  whenever  all  the  diagonal  elements 
In  a  reduced  matrix  are  zero,  l.e.  there  exists  no  permutation  N  such 
that  NAN*  has  an  L  D  L*  decomposition. 


2.9  Symmetric  Case  :  Failure  of  L  D  L*  In  Finite  Precision  Arithmetic 

The  L  D  L*  decomposition  in  finite  precision  can  be  unstable  if 
the  diagonal  elements  are  too  small.  The  L  D  L*  decomposition  on  Che 

R  fl 


matrix  A  « 


will  produce  the  same  Incorrect  solution  on  the 


i  ll 


computer  as  we  saw  In  §2,3  for  Gaussian  elimination  without  Interchanges. 
However,  here  there  exists  no  permutation  matrix  N  such  chat  NAN* 
has  a  stable  L  D  I,*  decomposition.  Thus  the  L  D  L*  decomposition 
falls  for  symmetric  indefinite  matrices. 


2.10  Symmetric  Case  :  Present  Situation 

If  we  Ignore  the  symmetry  of  A  and  apply  elimination  with  complete 
or  partial  pivoting  to  A  ,  then  in  general  A  will  no  longer  be  sym¬ 
metric  after  the  first  step  of  the  elimination.  We  then  need  n^ 

1  3 

Storage  positions  in  the  computer  and  we  must  perform  ^  —  n  multl- 
y  n’  additions.  But  we  also  have  stability  for  the 


plications  and 
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L  U  decomposition.  This  procedure  is  presently  recommended  for  Che 
solution  of  symmetric  Indefinite  systems  of  llnesr  equations  (Fox; 
p.  80,  185),  and  thus  the  symmetry  of  A  is  of  no  advantage. 

2.11  Symmetric  Case  :  Our  Problem 

A  is  symmetric  but  indefinite,  we  would  like  Co  find  an 

algorithm  which  gives  a  stable  decomposition  when  applied  to  N  A 

where  M  is  a  suitable  permutation  matrix,  but  which  also  takes 

advantage  of  the  symmetry  in  order  to  require  only  ~  y  storage 

1  3 

positions  in  the  computer  and  to  require  only  n  multiplications 

and  ~  4  n*  additions. 

6 

Our  algorithm  will  fulfill  all  the  above  requirements  with  Che 
exception  that  we  will  need  between  ~ n*  end  additions. 


Chapter  3  :  Historical  Survey 


3 . 1  IntroJuctlon 

Various  methods  have  been  proposed  tor  syaanetr.fc  Indefinite 
mairlces.  Most  of  these  methods  have  been  unstable,  while  the  stable 
methods  have  required  operation  counts  of  at  least  n^/3,  where  an 
operation  Is  defined  to  be  a  multiplication  followed  by  an  addition. 

Let  us  look  at  some  of  these  methods. 

^ 2  Direct  Methods 

Usually  direct  methods  tor  symmetric  Indefinite  systems  are  based 
on  the  symmetric  form  of  Gaussian  elimination,  L  Q  L  ,  which  is  unstable 
in  the  absence  of  pivoting. 

The  L  D  method  and  Its  variant,  the  method  of  congruent  trans¬ 
formations  (In  which  we  use  the  largest  diagonal  element  as  the  pivot 

1  »  Is 

at  each  step  (§2.7))  require  storage  locations  and  ^  " 

operations.  But  both  methods  are  unstable  (12.7-2.9).  These  methods 

are  of  value  only  if  It  is  known  in  advance  that  no  element  of  D  will 

vanish  or  be  small. 

The  Grout  factorization  (Hildebrand,  pp.  429-43S;  Householder, 
pp.  82-83)  is  also  a  modification  of  L  D  ,  symmecrlc  Gaussian  all- 

1  2  1  3 

mlnatlon,  and  thus  requires  n  storage  locations  and  ** 

operations,  but  It  Is  also  unstable. 

These  variants  of  L  D  are  unstable.  Let  ua  consider  some 
direct  methods  that  are  not  baaed  on  L  D  L^  . 


15 


The  escalator  method  (Householder,  pp.  78-79)  uses  Che  known 
solution  oC  a  subsystem  as  a  step  in  solving  the  complete  system. 

The  problem  lies  in  finding  the  solution  of  the  subsystem. 

Some  headway  in  this  problem  was  made  by  Parlett  and  Reid  (1969). 

They  reduce  a  symnetrlc  matrix  Co  tridiagonal  form  by  scabillzed  ele¬ 
mentary  congruences  and  solve  the  cridlagonal  system  by  Gaussian  eli¬ 
mination  with  partial  pivoting.  They  require  storage  loca¬ 

tions  and  ~  n^  operations,  and  Che  method  Is  stable. 

In  October,  1965,  U.  Kahan  (in  correspondence  with  R.  De  Meersmana 
and  L.  Schotamans)  proposed  a  method  for  solving  symmetric  systems  based 
on  Lagrange's  theorem  on  the  reduction  of  quadratic  forms  to  diagonal 
forms  (§4.3  -  4.4).  Kahan  proposed  the  generalization  of  the  idea  of 
a  pivot  to  Include  2^2  submaCr Ices  (§4.5). 

Then  Schotsmans  (1965)  prepared  an  algorithm  In  which  one  searches 

all  Che  principal  2^2  submacrlces  for  Che  one  with  largest  determinant. 

1  j  la 

This  algorithm  requires  "" 'J  **  storage,  but  between  ~ 

~ -I-  n’  operations  (§5.4). 

3 . 3  Indirect  Methods 

The  Seidel  iterative  method  (Householder,  pp.  48-51,  81)  and  the 

i 

method  of  relaxation  (Householder,  pp.  48-51,  81)  require  ~  n^  opera¬ 
tions  for  each  cycle.  The  number  of  cycles  required  depends  on  the 
matrix,  Che  starting  values,  and  Che  needed  accuracy >  Usually  the  num¬ 
ber  of  cycles  exceeds  n  ,  so  at  least  n^  operations  are  required. 
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The  method  of  steepest  descent  (Householder,  pp.  47-51,  82) 
requires  ^  2  n'  operations  at  each  step.  Again,  usually  at  least 
n  steps  are  required.  So  the  number  of  operations  Is  at  least 
^  2  n^  . 

The  congruent  gradient  method,  also  called  Che  Stlefel-Hescenes 
method,  is  a  finite  iterative  method  designed  for  positive  definite 
matrices,  but  it  can  be  used  for  symmetric  matrices  (Fox,  pp.  208-213; 
Householder,  pp.  73-78,  82).  It  requires  ~  2  n*  operations  at  each 
of  n  steps.  Once  again  2  n^  operations  are  required.  See  Reid 
(1967)  for  useful  observations  on  this  method. 


r 


Chapter  4  :  Diagonal  Pivoting 


4 . 1  Preaerving  Symeetry 

In  order  to  l.ave  a  direct  soethod  for  ayimaetrlc  matrices  which 
will  preserve  aymmetry,  we  can  perform  only  congruences  on  the 
matrix  A  ,  l.e.  If  we  prenultiply  A  by  a  non-singular  matrix  X  , 
then  we  must  also  poatmultlply  by  , 


4 . 2  The  L  D  L  Decompoaltion 

Let  ua  consider  Che  L  U  method  In  greater  detail.  We  con¬ 
vert  A  to  diagonal  form  by  congruences.  Let  us  consider  the  first 
step  of  the  decomposition: 

r  * 

a  C^l  a  0  ^ 

Let  A  - - .  If  a  /  0  ,  then  A  ■  L - - —  L  , 

C  B  0  B-CC  /a 


where  L  - 


C/a  I 


and  1  ,  la  the  identity  matrix  of  order  n  -  1 

n-1 


The  variant  of  L  D  L  called  Che  method  of  congruent  cranaforma- 
Clona  (Westlake,  p.  21;  De  Meersmans  and  Schotsmans)  uses  the  largest 
diagonal  element  as  pivot.  This  Is  equivalent  to  the  L  D  decom¬ 
position  of  N  A  ,  where  N  Is  a  permutation  matrix. 

As  we  saw  In  §2.8  both  methods  are  unstable  for  symmetric  (inde¬ 
finite)  systems. 

This  Instability  results  from  our  being  unable  to  bring  an  off- 
dlagonal  element  into  the  pivotal  position.  Since  we  are  using  only 
congruences,  we  can  bring  only  a  diagonal  element  Into  the  pivotal 
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position.  Wl^en  some  off-diagonal  element  ,  j  >  1  ,  Is  very  large, 

we  can  bring  Aj^  Into  the  (2,1)  position  by  congruences,  but  never 
into  Che  (1,1)  position.  Thus  we  cannot  take  advantage  of  this 
valuable  information. 


4 . 3  Orthogonal  Reduction  to  Diagonal  Form 

Any  real  quadratic  form  A  x  of  rank  r  can  be  reduced  by  an 
orthogonal  transformation  to  a  diagonal  form 

+  . .  +  , 

where  . X^  are  the  non-zero  eigenvalues  of  A  (Mirsky,  pp.  362- 

363). 


If  A  Is  an  n^n  symmetric  matrix  with  det  A  j*  0  ,  then  the 
above  means  that 

A  =  0  A  0*^  . 

where  A  »  diag  X,,..,,X  ,  the  X  are  the  eigenvalues  of  A  ,  and 

In  1 

0  la  an  orthogonal  matrix  whose  1^^  column  Is  an  eigenvector  corres¬ 
ponding  to  Xj  . 

However,  this  0  A  0^  decomposition  Involves  more  work  than 
Gaussian  elimination  and  requires  the  use  of  Irrational  operations. 

For  a  finite-precision  algorithm  we  would  prefer  a  reduction 
Involving  only  rational  operations. 


4.4  Lagrange ' s  Method  of  Reduction 

In  1759  Lagrange  devised  a  method  for  reducing  a  real  quadratic 
form  of  rank  r  by  a  real  non-singular  linear  transformation  to  a 
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where  are  all  non-zero  (Mlrsky,  pp.  368-374),  and  the 

number  of  positive  (and  negative)  squares  is  invariant. 

This  method  corresponds  to  the  L  D  decomposition  of  a  sym¬ 
metric  matrix  A  when  the  L  D  decomposition  exists. 

Suppose  A, ,»...«  A  “  0  ,  while  det  A  0  .  Then  the 
11  nn 

L  D  decomposition  for  A  does  not  exist.  In  this  case,  some 

A  ft  0  for  r  ft  s  since  det  A  i*  0  . 
tr  8 

► 

Let  us  assume  A^^^^  -  0  -  Ajj  but  >  where  A  -  A  . 

4>(x,,...,x  )  -  A  X  Is  a  quadratic  form  in  . . x  where 

in  in 

X*'  -  lXj^,...,X^]  . 

In  this  case  Lagrange  proposed  the  following  transformation: 


(4.4.1) 


yi  +  >2 


^1  '  ^2*  *3 


This  maps  2  Aj^2  *2  ^  ^12  ^^1  ~  ^2^  ‘ 

transformed  into  a  quadratic  form  <p  in  where  the  coeffi¬ 

cients  of  y^  and  y^  are  non-zero.  Then  we  can  proceed  with  the 
decomposition  (Mlraky,  p.  371-2;  Gantmacher,  p.  199). 


Let  us  consider  the  above  transformation  in  matrix  form.  Then 


T  y  "  ,  where 


(4.4.2)  T 


and  ^jj_2  Identity  matrix  of  order 


n  -  2  . 
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Thus  A  X  -  (T*"  A  T)  y  ,  and  AT  .b  a  symmeCrlc  matrix 

with  (T^  A  a  T)22  f*  0  • 


Hence  we  have  avoided  the  problem  of  zeros  on  the  diagonal  of  A 

■  1  1* 

by  use  of  the  2x2  matrix 

L  1  “1. 

This  procedure  Is  also  applicable  to  complex  quadratic  forms. 


A. 5  Kahan's  Proposal 

In  1965  W.  Kahan  (In  correspondence  with  R.  De  Meersmans  and 
L.  Schotamans)  proposed  that  Lagrange's  method  could  be  made  Che  basis 
of  a  sf’ble  method  which  preserved  symmetry. 

Kahan  adapted  Lagrange's  method  to  finite  precision  by  observing 


chat  Che  use  of 


in  (A. A. 2)  corresponds  to  the  use  of  a  2^2 


submatrlx  as  a  pivot  in  a  decomposition  by  linear  transformations  and 
that  a  2x2  could  be  chosen  if  the  diagonal  elements  were  zero  or  very 
small  (§§2.8-2. 9). 


Suppose  we  used  a  2x2  submatrlx  P  as  a  pivot.  Let  us  look  at 


such  a  decomposition.  Let  A  > 


> 

_c 

B 

,  where  A  -  A  ,  det  A  0  , 


A  is  nxn  ,  P  is  2x2  ,  C  is  (n-2)x2  ,  snd  B  is  (n-2)x(n-2).  Then 


A  -  L, 


P 

0 

L^  ,  where  Lj^  ■ 

0 

B-CP~^ 

CP"^ 

I  . 

J 

n-^ 

and  I,  is 
k 


the  identity  matrix  of  order  k  . 

How  do  we  choose  whether  to  use  a  1x1  or  a  2x2  pivot?  Can  both 
1x1  and  2x2  pivots  be  bac? 
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4.6  Pivotal  Strategy 

Kahan  considered  two  pivotal  strategies.  In  the  first,  one 

searches  the  entire  oiatrlx  for  m  ■  sax  {A 

l.J.k 


2 

11 


|A 


Jj 


If  "“c  *  Interchange  rows  and  coluans  1  and  1  and  use 

aa  a  1x1  pivot.  If  j I  •  then  Interchange  rows  and 


columns  1  with  j  and  2  with  k 


and  use 


as  a  2x2  pivot. 


Since  we  search  the  entire  matrix,  this  Is  called  a  complete 

pivoting  strategy.  In  analogy  with  complete  pivoting  for  Caussian 

Elimination.  However,  the  searching  here  requires  between  ^  ~  n^ 

and  n’  multiplications  to  find  for  all  steps  (depending  on 

the  number  of  ixl  and  2x2  pivots  used).  The  decomposition  Itself 

requires  ~  ^  n’  multiplications.  Thus  this  strategy  would  require 
Is  Is 

between  ~  -j  n  and  ~  ”  multiplications  to  solve  A  x  ■  b 

(§3.4),  which  Is  more  than  for  Gaussian  elimination.  Hence  Kahan 
rejected  this  strategy. 

Kahan  considered  a  second  pivotal  strategy  In  which  we  scan  only 

the  first  column  and  the  main  diagonal;  this  Is  called  a  partial  pivoting 

strategy,  in  analogy  with  partial  pivoting  for  Gaussian  elimination. 

1  2  1  2 

The  searching  here  only  requires  between  ~  ^  ~  *2  "  multi¬ 

plications  . 

We  take  m  “  max  ^A*^,  |Aj^j^  A^^  -  A^]}  .  However,  this  partial 

1  >  j 


pivoting  strategy  is  unstable. 
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Ut  A 


e 

e 


c  c 

^  l' 
1  ¥ 


,  where  0  <  e  «  1  . 


Then  n 


1^11  ^22 


'21' 


Thus 


1 

^  e 
e  ^ 


would  be  used  as 


/jN  2  1  8  2 

a  2^2  pivot,  and  the  reduced  naCrlx  a'’  T  e  ^  T  ~  3  ‘ 

If  c  Ih  small  enough,  then  in  finite  precision  arithmetic  the 


2  1  8  2  2  1  8 

operation  yields  “  3  ^  ^  3  '  cause  highly 


inaccurate  solutions,  as  in  §2.3. 

So  this  partial  pivoting  strategy  is  unstable.  For  these  reasons 
Kahan  rejected  this  method  for  use  on  symmetric  systems. 


4 . 7  Parlett's  Observation 

In  1967  B.  Farlett  observed  that  the  examples  for  which  the  partial 
pivoting  strategy  was  unstable  were  also  unequilibcated .  A  symmetric 
matrix  A  is  equilibrated  if  max  |a. .|  ■  1  for  each  row  index  1 

J  ^ 

(§7.1).  Parlett  conjectured  that  the  partial  diagonal  pivoting  strategy 
would  be  stable  when  applied  to  equilibrated  matrices. 
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Chapter  5  :  The  Decomposition  for  PtaKonal  Pivoting 


5.1 


Definitions 


Let  A  be  an  nxn  ayimetrlc  non-alngulai  matrix.  We  want  to 
reduce  A  to  the  "diagonal"  form  M  D  by  congruences,  where  D 
la  a  block  diagonal  matrix,  each  block  being  of  order  1  or  2,  end  M 
la  unit  lower  triangular  with  i  *  ®  *^1+1  i  ^  ' 

Let  Uq  “  MX  Ia^jI  ,  -  max  |a^^1  ,  and  v  -  |a^j^  -  A^^^j  , 


5.2 


l.j 

The  Decompoaltlon 


Let  A 


P 

c^‘ 

_c 

If  ?  ^  exists,  then  A  L. 


Is 

J  (n-J), 

2. 

P 

0 

p 

B-C  P”^  C*^^ 

L^  ,  where 


L  • 

1 


C  P 


-1 


n-jj 


,  and  I, ,  I  ,  are  the  Identity  matrices  of  order 
J  n-J 


,-l 


j  and  n-J ,  respectively.  Any  element  of  CP  will  be  called  a 
multiplier. 


5.3 


The  1x1  Pivot 


Suppose  P  is  of  order  1.  (We  shall  not  make  a  distinction 
between  a  matrix  P  of  order  1  and  Its  element,  which  we  shall  also 


call  P  .) 


2A 


Let  us  assume  that  we  have  already  interchanged  rows  and  columns 
so  that  |p|  <•  ,  l.e.  P  Is  the  maximum  diagonal  element. 

If  P“^  exists  (l.e.  ^  0),  the.,  let  -  B  -  C  p“^  . 

Then  (CP  )^  -  A^^j^/P  and  A^j  -  ^1+1, J-H  ~  ^  ^1  ^j+1,1  ‘ 

Since  |p|  -  y  and  U-  -  max  |a. . |  ,  we  have  the  following: 

l,j  ^ 

Lemma  1:  If  P  la  of  order  1  and  |P|  ■  /  0  ,  then 

(1)  max  1(C  1  Mq/Mj  , 

(il)  max  1a^j“^^|  <  (1  +  . 

j 

Thus  a  1*1  pivot  P  Is  useful  Iff  |pj  “  hj  la  large  relative 
to  y^  ,  l.e.  If  Is  bounded  away  from  zero. 


5.4  The  2x2  Pivot 

Suppose  P  Is  of  order  2  and  F  ^  exists  (l.e.  v  0) . 

Here  the  (k-1)®^  row  of  C  P  ^  is: 

^\l’  ^2^  ^  “  All  ^22  ■ 

Ut  A<"-»  .  B  -  C  p-‘  C'  .  Then  -  (C  p-h^i  Cjj  - 

(c  . 

Since  V  “  1^11  ^22  ”  ^21^  ’  ^^1®  ^^2^  —  ^0  *  a***^ 

|Aii1,  Ia.221  1  Ml  ,  we  have  the  following: 
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Leana  2: 


(1)  max 

i.J 

(11)  max 
i.j 


If  F  la  of  order  2  and  |det  P|  •  V  0 
|(C  <  Mq  %  +  V^)/^  , 

<  (1  +  2Uq  (Uq  +  Uj^)/v]  Uq  . 


then 


Thus  a  2^2  pivot  is  useful  Iff  we  can  bound  V  away  from  zero. 
In  particular  from  §5.3,  we  need  to  have  v  bounded  away  from  zero 
ifhenever  la  near  zero. 


(Note  that  Che  use  of  the  standard  norm  bound  would  give  too  crude 
an  analysis  for  Lemmas  1  and  2.) 


5 . 5  Bounding  u 

We  can  easily  bound  v  from  above,  since  v  ■  -  A^^l  < 

|A2i|^  +  IAjj^I  Ia^jI  1  Uq  +  Uj  .  Thus  we  have: 

Leama_3:  |det  P  |  >•  V  5.  Mq  +  • 


This  upper  bound  Is  sharp  for 


But,  as  we  saw  In  §5.4 


we  need  a  lower  bound  on  v  which  bounds 
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V  away  £rou  zero  when  VIj^/Uq  Is  soiall.  Clearly,  such  a  lower  bound 
does  not  exist  without  Interchanges. 


We  shall  exhibit  three  different  pivotal  strategies  In  §6.1,  §6.4, 
and  §8.1,  which  provide  us  with  the  necessary  Interchanges  so  that  we 
have  a  2^2  pivot  P  with: 

Idet  P|  »  V  1  p*  -  p*  .  (§6.2,  §6.5,  §8.2). 

Assuming  this  lower  bound,  we  have: 

Lemma  4;  If  |det  P|  -v>u^-Pj>0  ,  then 

(1)  max  |(C  for  1^-1. 2, 

1  ,Ic 

(11)  max  £  Wq  [1  +  2Uq/ (Pq  -  Pj^)]  . 

The  lower  bound  on  V  >  0  Is  sharp  for 


Thus  we  have  a  good  bound  on  element  grouch  in  the  reduced  matrix, 
since  if  Vj^/Pq  is  small  then  1  -  Uj^/Pq  +  1  •  We  shall  see  in 
Chapter  10-12  that  stability  follows  from  this. 
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5.6 


The  Reduced  Matrices 


Define  A 


(n) 


/n) 

0 


Ul  .  V 


(n)  _ 


Let  A 


(k) 


be  Che  reduced  matrix  of  order  k 


l.j  J 


-  max  1a;'‘>|.  and  v<‘^> 
i  “ 


Let 

iAk)  (k) 
'*11  *22 


-  A 


(k): 

21 


I  . 


All  considerations  In  SSS.2-S.5  hold  for  A 


(k) 


5.7  Criterion  for  Choosing  a  1X1  or  2><2  Pivot 

Ue  must  find  a  proper  criterion  for  deciding  whether  we  shall  use 
a  1^1  or  2x2  pivot. 

In  Chapter  10  we  shall  show  that  the  elements  of  the  error  matrix 
are  bounded  In  proportion  to  the  elements  In  the  reduced  matrices.  For 
stability  we  must  ensure  that  the  elements  in  the  reduced  matrices  do 
not  become  Coo  large. 

If  we  made  our  criterion  to  be  the  minimization  of  the  number  of 
mvtlt Ip Ilea clone  (additions),  then  we  would  want  a  1^1  (2^2)  pivot  at 
each  step.  But  this  would  be  unstable. 

Instead,  let  us  aim  to  minimize  the  element  growth  chat  can  take 
place  In  the  transformation  from  one  reduced  matrix  to  the  next.  For 
further  remarks  see  S12.6. 


Let 

pivot  for 


be  the  growth  factor  permitted  by  choosing  a 
,  where  j-1  or  2. 


If  the  hypothesis  of  Lemma  4  holds  (l.e.  ^ 

(k) 

for  all  ,  Chen  by  l-ensias  1  and  4: 
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,  -  1  +  2/(1  -  y} 


has  a  good  bound  If  is  not  coo  small;  while 

Ck^  (k^  (k^ 

Fj  '  has  a  good  bound  If  /Uq  la  not  too  large.  Thus  we  are 

led  to  Che  following: 

Definition:  For  0  <  a  <  1  ,  let  be  the  following  strategy: 

(k) 

for  each  reduced  matrix  A  ,  choose  a  1x1  pivot  iff 

(^\  /’Irl 

^1  ^^0  pivot  otherwise). 


we  have 


<  1  +  1/a 


^2  1  1  +  2/(1  -  a)  for 


all  A' 


But  at  any  stage  the  choice  of  a  2x2  pivot  carries  us  further 

cowards  Che  complete  reduction  than  does  the  choice  of  a  1x1  pivot. 

Since  Che  growth  factors  from  reduced  matrix  to  reduced  matrix  are 

multiplicative,  it  is  natural  to  compare  the  square  of  the  growth  fac- 

(k)  (k) 

tor  Fj^  permitted  by  choosing  a  1x1  pivot  for  A  with  the  growth 

(k) 

factor  '  for  a  2x2  pivot. 

Thus  the  problem  Is  to  find  min  max  {(1  +  1/a)*,  1  +  2/ (1  -  a))  . 

0<a<l 

Theorem:  min  max  {(1  +  1/a)*  ,  1  +  2/(1  -  a))  ■  (9  +  »'TT)/8 
0<a<l 

and  is  achieved  by  a  -  Oq  -  (1  +  i/lT)/8  . 

Proof :  The  equation  (1  +  1/a)*  -  1  +  2/ (1  -  a)  reduces  to  a 
quadratic  with  roots  (1  ±  i'^l7)/8  .  Since  the  left  side  of  the  equation 
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is  monotone  decreasing,  Che  right  aide  is  monotone  increasing,  and 
a  >  0  ,  the  minimum  is  given  by  a  «  -  (1  +  /r7)/8  •  q.e.d. 


We  immediately  obtain  the  following  bounds  on  multipliers  and  on 
elements  in  Che  reduced  matrices  under  strategy  .  Let  m  be  any 
multiplier. 


(Ic) 

Corollary  2:  Under  strategy  ,  if  for  all  A  , 


Corollary  3;  Under  strategy  with  a  -  ,  if 

for  all  ,  then  for  1  1  i  1  n  : 

i  ^^oU9  +  <  Uq  (2.57)""^  . 

In  Chapters  11-12  we  will  give  a  much  better  bound  on  the  elements 
in  the  reduced  matrices. 

The  strategy  allows  us  to  proceed  in  the  following  order; 

(1)  calculate  and  ; 

(k)  (k) 

(2)  if  ^0  ®  pivot; 

(3)  otherwise  we  find  a  2x2  by  some  strategy. 

Thus  we  search  for  a  2x2  for  iff  oig  • 
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Chapter  6  :  The  Complete  and  Partial  Pivoting  Strategies 

6.1  Complete  Pivoting 

we  saw  In  Chapter  4,  the  partial  pivoting  strategy  can  be 
unstable  when  used  on  unequlllbrated  matrlcis.  The  trouble  lies  in 
the  fact  Chat  we  do  not  have  a  lower  bound  on  the  2^2  principal  minors 
which  bounds  them  away  from  zero  when  the  diagonal  elements  are  small 
(§5.5). 

Let  us  therefore  consider  a  complete  pivoting  strategy  ("complete" 
in  the  sense  that  we  search  over  all  the  principal  2^2  minors,  cf.  §4.6). 

By  Interchanging  rows  and  the  corresponding  columns  It  is  possible 
to  bring  any  diagonal  element  into  the  (1,1)  position  or  any  principal 
2x2  submatrix  into  the  leading  2x2  position. 

Let  ^  Ajj  -  Ajj!  . 

^  »  J 

The  complete  strategy  involves: 

(1)  finding  U,  -  max  |a  J  ,  u-  -  max  |a  j  ; 

1  l,j  ^ 

(2)  choosing  a  1x1  or  2x2  pivot  according  to  (§5.7); 

(3)  for  a  1x1,  interchanging  so  that  jp]  -  (§5.3); 

(4)  for  a  2x2,  finding  ,  and  interchanging  so  that  |det  P|  • 

(§5.4). 

This  would  be  repeated  for  each  reduced  matrix. 
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6 . 2  Bound  InR  \)^ 

However,  the  result  of  all  this  work  Is  that  we  do  obtain  a  lower 
bound  for  in  terras  of  Uq  and  . 

Theorem  1:  Pq  -  1  <  p*  +  p[. 

Proof:  The  upper  bound  follows  from  Lemma  3  of  §5.3.  Let 


"  l\sl 


Then  ^ 

^  »  J 


A* J  >  |a  a  -  A^  I 

li '  —  •  rr  88  rs ' 


Un  -  A 


rr  88 


i  Pq  -  pj  .  since  Pq  -  A*^  and 


Pi  -  max  lA^J 


q.e.d. 


Thus  Lemma  4  in  §5.5  holds  for  when  p^^  <  p^  .  According  to 

(kl 

,  we  choose  a  1x1  pivot  for  the  reduced  matrix  A  if 

pj*"^  >  Qq  pJ*"^  .  where  -  (1  +  .^)/8. 

In  Chapters  10-12  we  shall  see  that  stability  of  this  complete 
pivoting  strategy  follows. 


6.3  Operation  Count  for  Complete  Pivoting 

Unfortunately,  the  operation  count  here  is  much  larger  than  we 
desire . 

(k) 

The  calculation  of  requires  k(k-l)  multiplications  and 

additions.  Let  p  be  the  number  of  ixl  pivots  used  (so  q  »  (n  -  p)/2 

pivots  of  order  2  are  used).  Let  denote  summation  over  those 

(kl 

indices  k  ,  1  £  k  ^  n  ,  for  which  a'  ' 


uses  a  pivot  of  order  j  . 
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The  searching  for  all  the  requires: 


(6.3.1)  k(k  -  1)  +  k(k  -  1)  multiplications  and  additions. 

n  I  . 

Now  (6. 3.1)  £  J  k(k  -  1)  n(n  +  1)  (n  -  4)  ~  rr  n*  with 
r-1 

1 

equality  iff  p  -  n  ,  while  (6.3.1)  ^  ^  2J  (2j  -  1)  ■  Yo  2)  (2n  “  1) 

j-1 


-  i  _3 


n*  with  equality  Iff  p  -  0  (l.e.  q  ■  n/2). 


t'ruto  Chapter  9  we  see  that  the  rest  of  the  work  for  solving  A  x  ••  b 
would  require  ~  ^  n^  multiplications,  and  between  ~  ^  n^  and 

~  additions. 

Thus  the  complete  diagonal  pivoting  strategy  requires  between 
~  •j  and  ~ n’  multiplications,  and  between  and  ~ -j  n* 

additions.  (To  be  exact,  n^  ^  multiplications  and 
~  n*  ^  additions  are  required,  where  p  is  the  number  of  l^l 


pivots  used.) 


Examples:  If  A  is  n’<n  and  positive  definite,  then  p  ~  n  . 


If  A 


is  nxn  with  n  even,  then  p  “  0 


0 


■1 
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6.4  Partial  Pivoting  for  Equilibrated  Matrices 

As  we  saw  In  §§6.2-6. 3  the  complete  pivoting  strategy  is  stable, 
but  it  Involves  more  work  that  we  are  willing  to  perform,  while  In 
Chapter  3  we  saw  that  a  partial  pivoting  strategy  Is  unstable  for 
unequlllbrated  matrices. 

We  shall  now  show  that  a  partial  pivoting  strategy  is  stable 
when  applied  to  equilibrated  matrices  (with  equilibrated  reduced 
matrices) . 

Let  A  be  equilibrated,  l.e.  let  max  ■■  Wq  for  every  1, 

where  >  0  (usually  we  normalize  by  taking  Uq  *  1  ) • 


Let 


V  ■■  max 

P  I 


Kii 


The  partial  strategy  Involves: 

(1)  equilibrating  A  (thus  we  know  ); 

(2)  finding  and  choosing  a  ixl  or  2x2  according  to  (§5.7); 

(3)  for  a  Ixl,  lnterchai\ging  so  that  jp]  ■  p^^  (§5.3); 

(4)  for  a  2x2,  finding  ,  and  Interchanging  so  that  |det  p]  « 
(§5.4). 


6 . 5  Bounding 

Now  let  us  find  a  lower  bound  for 
Theorem  2;  If  A  Is  equilibrated 

then  -  Uj  1  Vp  <  p^  +  pj  . 


V  . 
p 

(max 

J 


Pg  for  every  i  ), 


Proof ;  The  upper  bound  follows  from  Lemma  3  of  §5.5  By 


equilibration,  either  (i)  there  exists  integer 

k  _>  2  with  ^0  "  ^1 

Vp  i  -  P*  -  0  .  If  (ii),  then  Vp  >  !Ajj^  -  A^J  - 

^0  -  \l 

Thus,  l^oinia  4  in  §5.5  holds  for  if  .  According  to 

(It) 

,  we  choose  a  1x1  pivot  for  the  equilibrated  reduced  matrix  A 

of  order  k  iff  ^  ,  where  <Xq  ■  (1  +  /r7)/8  . 

Also,  stability  of  the  partial  pivoting  strategy  for  equilibrated 
matrices  follows  from  Chapters  10-12. 

For  A  «  A^”^  ,  only  2(n  -  1)  multiplications  and  additions  are 

required  to  calculate  .  Thus  the  calculation  of  for  all 

P  P 

1  2 

k  requires  between  n  and  n(n  -  1)  multiplications,  compared 

with  between  ~  T the  in  (6.3.1). 

0  0  c 

6 . 6  Criticism  of  the  Partial  Pivoting  Strategy 

The  drawback  to  this  method  and  the  criterion  which  we  have  found 
for  the  pivoting  strategy  is  that  the  matrix  mus::  be  equilibrated  at 
Che  beginning  and  then  each  reduced  matrix  sh.ould  be  equilibrated. 

But  an  algorithm  for  equilibrating  symmetric  niatrlces  has  never  been 
exhibited,  l.e.,  for  an  arbitrary  matrix  A  ,  vie  seek  diagonal  matrices 
Dj^  ,  such  that  Dj^  A  is  equilibrated,  .nd  if  A  is  symmetric, 

we  need  »  D2  in  order  to  preserve  symmetry . 


35 


We  resolve  the  above  predicament  by  two  fundamentally  different 
approaches.  In  Chapter  7  we  exhibit  an  algorithm  which  can  equili¬ 
brate  any  symietrlc  matrix  in  a  very  simple  way.  In  §§6.1-6. 3  we 
showed  Chat  complete  diagonal  pivoting  avoids  the  problem  of  equili¬ 
bration  and  is  stable,  although  the  number  of  multiplications  and 
additions  required  is  more  than  we  desire.  But  in  Chapter  8  we  shall 
exhibit  a  new  version  of  diagonal  pivoting  which  is  applicable  to 
unequilibrated  matrices;  we  call  this  unequilibrated  diagonal  pivoting. 
Ihia  method  will  show  that  equilibration  (in  Wilkinson's  sense)  is 
unnecessary  for  this  strategy  and  this  is  the  algorithm  that  we 


recommend. 
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Chapter  7  ;  Equilibration  of  Sympetric  Matrlcea 

7 . 1  Introduction 

Wilkinson  (1961)  recotoDends  that  a  matrix  be  equilibrated  before 
applying  any  algorithm  for  solving  a  syatem  of  linear  equations.  A 
matrix  A  is  said  to  be  equilibrated  if  all  its  rows  and  columns  have 
the  same  length  in  some  norm.  Wilkinson's  rounding  error  analysis  for 
Gaussian  elimination  (Wilkinson,  1961)  gives  the  most  effective  results 
when  the  matrix  is  equilibrated,  since  a  small  perturbation  of  one  row 
(or  column)  Is  then  of  the  same  magnitude  as  that  of  any  other  row 
(or  column). 

7 . 2  Equ^ libration  of  General  Matrices 

In  finite  precision  we  modify  the  definition  of  equilibration. 

(In  this  chapter  we  shall  confine  ourselves  to  the  norm 

lxl„  -  .) 

A  matrix  A  is  row  equilibrated  if,  for  each  row  index  i  , 

0  <  max  |6. ,1  ^  Wq  .  where  3  is  the  number  base  of  the  floating 

point  system.  A  matrix  A  is  column  equilibrated  if,  for  each  column 

index  J  .  B  ^  Pf,  5.  max  |a  .  |  <  p^  . 

l<J<n  J 

A  matrix  A  is  equilibrated  if  it  is  both  row  and  column  equili¬ 
brated  . 

Usually  we  normalize  by  choosing  Pq  -  1  .  We  shall  assume 


1  in  this  chapter. 


The  use  of  6  permits  a  matrix  to  be  equilibrated  by  changes  of 
exponent  only. 

In  order  to  row  (column)-equlllbrate  A  we  seek  a  diagonal  matrix 
(D^)  such  that  A  (A  D2)  is  row  (column)-equlllbrated .  To 
equilibrate  A  we  seek  diagonal  matrices  that  A  Dj 

la  equilibrated.  Ikswever,  there  Is  no  unique  equilibrated  form  of  a 
matrix  for  this  norm  (Forsythe  and  Holer,  p.  45).  The  various  equili¬ 
brated  forms  may  differ  greatly  In  their  desirability  for  use  in 
Gausalan  elimination,  since  the  various  equilibrations  cause  different 
choices  of  the  pivots  (some  are  good  choices,  others  are  bad). 

I,  _  _  /II  ?  !  X .  I  )  the  equilibration  is  unique 

For  the  one-norm  (Ixl.  ■  ^  '  I'  ^ 

^  1-1 

but  the  convergence  of  the  algorithm  is  slow  (Slnkhorn  (1964),  (1967); 
Slnkhorn  and  Knopp). 


7 , 3  Difficulties  VII th  Symmetric  Matrices 

-2  -1 

If  A  ia  aymmt' <.x  Ic ,  then  we  allow  6  Uq  instead  of  6 
as  the  lower  bound  li  the  definitions  of  row  (and  column)-equllibrated . 
If  A  ia  syametrlc,  then  A  is  row  equilibrated  Iff  A  Is  column 
equilibrated  'iff  A  is  equilibrated. 

In  order  to  equilibrate  a  symmetric  matrix  and  still  preserve  the 
symmetry  we  seek  a  diagonal  matrix  D  such  that  DAD  is  equilibrate 


Let  A 


1  2 
2  1/2 


-  dlag  {1/2,  1}  ,  D^ 


diag  {/2/4,  /2  }  . 


dlag  {1,  1/2}  ,  and 
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Then  0^  A 


1/4  1 

1  1/2 


1  1 


,  D,  A  D  -  ,  and 

^  ^  11/8 


D3  A  D3  =  1^  ^ 


[1/8  1 


are  all  equilibrated. 


Criteria  for  choosing  pivots  are  based  on  their  size.  It  Is  an 
open  problem  whether  all  equilibrated  forms  of  a  symmetric  matrix  give 
satisfactory  pivots  as  Judged  by  the  usual  criteria. 

There  is  no  known  algorithm  for  the  symmetric  case. 


7 . 4  The  Obvious  Attempt 

Let  us  cc  slder  DAD.  Let  D  ■  dlag  {d,  , . . .  ,d  }  and  A  be 

1  n 

an  n'^n  symmetric  matrl-X.  Then  (DAD)^^  ■  d^^  A^^  .  Let  us  assume 

that  no  row  of  A  Is  all  zero. 

The  method  which  seems  the  more  obvious  is  to  equilibrated  one 
row  at  a  time:  Let  D^  *  diagfd^,! , . . . ,l)  . 


for  i ,J  > 


Then  (D  AD  )  -J  d.  A  for  J  1,  Jj*! 

l^dj^  Aj^j^  for  i,J-l 


So  choose  “  “tax  |a,  .  |  ,  max  |a  |)  . 

Now  is  symmetric  and  max  |(D.AD.)^  1  -  1  .  However, 

^  l<J<n  J 


when  we  try  ro  equilibrate  the  second  row  we  usually  destroy  the  equi¬ 


libration  of  the  first  row. 
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Proof :  By  Induction.  Let  “  1//  ^  ® 

by  hypotheela.  So  |d*  “  1  • 

Assume  that  d^,...,d^  1  ^  ^  have  been  chosen  so  that 

max  |d  d  B  .  |  «  1  for  1  ^  j  1~1  • 

Let  d  ^  =  max  {  ,  max  jd  B  |  }  .  By  hypothesis, 

^  1^<1  J 

d  ^  >  0  .  Then  max  |d  d  “  1  • 

j  1  3  13 

Hence  the  d^  exist  for  3.  ^  1  ;<  n  by  Induction.  Let 
D  -  dlag  {dj^,..,d^}  ,  Then  D  B  D  is  row-equilibrated,  q.e.d. 

7.6  The  Algorithm  for  Null  Rows  in  the  Lower  Triangle 

But  what  do  we  do  if  Aj^j^  *0  or  if,  for  some  i  ,  A^j  »  0 
for  1  1  3  1  i  ? 

Let  us  form  D  as  in  §7.5  with  the  exception  that  we  set 
d^  =  1  if  =  0  for  1  <  3  <  i  , 

Thenfor  all  l,j  :  |  (D  A  l*)jj  I  £  1  • 

DAD  falls  to  be  equilibrated  only  if  for  some  i  : 

(a)  A^j  =  0  for  1  ^  3  1.  1  ,  or 

(b)  max  |d  A  |  <  1  . 

3>1  J  ^ 

If  the  i  row  of  A  is  not  null,  then  the  maximum  in  (b)  is 

positive.  For  such  i  ,  define  e.  by  e  ^  =  max  |d  A  |  ;  for  all 

^  ^  3>1  ^ 

other  1  ,  let  e^  =  1  .  Let  E  “  dlag  {ej^,...,e^}  . 
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Theorem ;  Let  A  be  an  n^n  symmetric  matrix  with  no  null  row. 

Let  D  and  £  be  constructed  as  above.  Let  A  D  E  .  Then  A  is 
a  diagonal  matrix,  and  AAA  is  equilibrated. 

Proof ;  If  e^  1  then  the  maximal  magnitude  of  row  1  in  DAD 
is  raised  to  1  by  forming  E  D  A  D  E  ,  while  in  all  other  rows  the  i*"^ 
element  is  increased  in  magnitude  but  not  in  excess  of  1  . 

The  theorem  follows  from  57. 5.  q.e.d. 


7 . 7  Summary  of  Equilibration  of  Symmetric  Matrices 

If  A  is  an  n^n  symmetric  matrix  with  no  null  row,  we  can  find 

a  diagonal  matrix  D  (in  two  sweeps,  although  only  n  steps)  such 

that  DAD  is  equilibrated  (in  the  “-norm,  ■  max  |x. |  ). 

i  ^ 

We  can  express  D  by  the  following  algorithm:  For  i  =  l(l)n  : 


max  {  / 


max 

1S)11-1 


|«J  Ayl) 


if  A^j  ^  0  for 
some  j  ,  1  5.  J  £  ^ 


1^1  if  -  0  for  1  1  J  1  i  . 

Then,  for  1  -  l{l)n: 


if  A^j  #  0  for  some  j  ,  f  J  f.  f 
A 


max  |6  A  j  if  A  -  0  for  1  j  <  1 

V.  i+l^<n  ^  J  J 

Let  D  “  dlag  {dj^,...,d^}  .  Then  DAD  is  equilibrated. 

The  work  required  to  equilibrate  an  n>«n  symmetric  matrix  with  no 
all  zero  row  is: 

n  square  coots 

n(n+l)  multiplications 

(n-1)' 


iidul  tlons 
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In  practice  the  algorithm  can  be  expressed  In  a  very  simple 
manner  In  Algol  (see  Appendix  B)  and  can  be  performed  In  only  n  steps. 
(We  would  actually  set  f *0  If  -  0  for  1  ^  j  ^  1 

then  search  for  max  (6  A  |  Iff  f  -  0  ). 

l+l<J<n  J  ^ 
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d.+d. 


_  4  4 

This  will  give  6  *  ^  max  |6  ^ 


1  ^  1  for  every  1 


(This  equilibration  can  be  performed  very  rapidly  In  roach  ne 


language . ) 


Chapter  8  :  Unequlllbrated  Diagonal  PlvotlnR 


8.1  Maximal  Off-diagcnal  Element 

In  order  to  obtain  a  lower  bound  of  ~  uj  for  v  in  Theorem 

U  i  p 

1  In  §6.3,  we  needed  the  fact  that,  due  to  equilibration,  there  existed 
an  element  of  maxlinal  absolute  value  in  the  first  column.  However,  if 
Vij  Uq  )  then  there  exist  Integers  l,j  with  i  >  J  ,  such  that 

I  “  Mq  ■  We  need  only  bring  the  element  up  to  the  (2,1) 

position  and  then  we  will  have  a  2*2  pivot  with  a  maximal  off-diagonal 
element.  We  shall  call  this  variation  unequilibrated  diagonal  pivoting. 
Let  Wq  «  max  IAj^j  |  ■  >  where  r,3  are  the  least  such  inte- 

1 1  j 


gers.  Let  t*].  “ 


*li' 


Let  V.  =  A  A 

b  '  rr  ss 


A-  I  . 
rs ' 


This  strategy  involves: 

(1)  finding  and  the  least  integer  k  with  “  hj 

(2)  finding  ii^  and  the  least  integers  r,s  with  “ 

(3)  choosing  a  l^l  or  2^2  pivot  according  to  (§5.7); 


^0  ' 


(4)  for  a  I’^l,  interchanging  rows  and  columns  1  with  k  so  that 

|P|  -  (§5.3); 

(5)  for  a  2^2,  Interchanging  rows  and  columns  1  with  r  and  2  with  s 


so  that  jdet  P|  =  and  1^21 ! 


Wq  (§5.4). 


This  procedure  is  repeated  for  each  reduced  matrix. 


Note  that  calculating  requires  only  2  multiplications  Instead 

of  n(n-l)  for  v  and  2(n-l)  for  v  . 

c  p 

Clearly  from  the  definitions  of  •  V  ,  and  v  and  from  Lemma 

b  p  c 

3  o  f  y  5 . 5 ,  we  have  ; 

Lemma  1:  V.  <  v  <  V  <  y?  +  y*  . 

b  —  p  —  c  —  0  1 


8.2  Bounding 


Let  US  now  bound  from  below.  From  §8.1,  we  may  assume 

\  ”  ^\l  ^22  ”  ^21 1  1^21^  *  ^0  * 

Theorem  1;  If  -  y^  ,  then  ~  ^*1  1  -  ^0  ^1  * 


Proof :  The  upper  bound  follows  from  I.enima  1. 

Since  IA21I  -  Mq  .  |  -  y^  -  A^^  A22  > 


u 


2 

1  ' 


Here,  as  in  §6.2  and  §6.5,  symmetry  was  used  to  get  the  lower  bound 
on  .  By  symmetry,  if  y^  =  1A2j^|  then  |Aj^2l  “  1*0  ‘ 

If  A  were  not  symmetric,  then  y^  =  |A2^|  does  not  imply  that 

|Aj^21  ”  Vq  (f*'  fact  we  could  have  Aj^2  "  Thus  no  such  non-negative 

lower  bound  on  the  determinant  of  2x2  submatrlces  can  exist  for  non- 
syometrlc  matrices. 

Thus  Theorem  1  implies  that  Lemma  A  in  55.5  holds  for  v.  if 

b 

y^  <  pQ  .  According  to  (§5.7),  we  would  choose  a  1x1  pivot  for  the 

fk)  fkl 

reduced  matrix  A  of  order  k  iff  y^  ^  y^  ^  ,  where 

Uq  =  (1  +  .T7)/8. 


In  Chapters  10-12,  we  shall  see  that  stability  of  this  strategy 


(for  unequillbraced  matrices)  follows  from  the  above. 
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8 . 3  Cofliments  on  this  Strategy 

From  §8.1,  we  see  that  we  need  do  no  searching  over  the  2^2  prin¬ 
cipal  minors,  but  we  merely  choose  that  principal  minor  with  maximal 
olf-dtagonal  element.  In  Chapter  9  we  shall  see  that  this  searching 
for  the  requires  between  ~  ^  additions 

and  no  multiplications. 

We  used  the  terms  "complete"  and  "partial"  strategies  in  Chapter 
4  to  distinguish  between  searching  over  all  the  2x2  principal  minors 
and  over  the  2x2  principal  minors  with  off-diagonal  element  in  the 
first  column. 

In  analogy  with  Gaussian  elimination  with  complete  pivoting,  we 

would  like  to  call  our  strategy  in  §8.1  a  complete  pivoting  strategy 

(k) 

since  we  search  all  of  A  for  Its  maximal  element,  but  we  do  not 

wish  to  cause  confusion  with  the  use  of  the  word  "complete"  in  Chapter 

4  where  it  meant  searching  all  the  2x2  principal  minors. 

Now,  in  analogy  with  Gaussian  elimination  with  partial  pivoting, 

we  ask  if  there  could  be  a  partial  strategy  where  we  search  only  the 
(k) 

first  column  of  A  for  its  maximal  element. 

Such  a  partial  strategy  requires  at  most  only  ^  n  (n-i.)  addi¬ 
tions  to  calculate  the  maximal  element  in  the  first  column  of  all  the 

reduced  matrices.  If  such  a  partial  strategy  were  stable,  then  this 

1  3  1  • 
strategy  would  require  only  ""  t  n  multiplications  and  ~  x 

t>  o 

additions  to  solve  Ax-b  ,a=A^  ,  detAj^O  . 
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However,  any  such  partial  strategy  la  unstable  for  unequlllbrated 
matrices,  as  the  following  example  shows: 


a  e  £  £ 

e  'a  ®  ^  ^ 

1 

t  e  1  2'  “  ^ 


where  0  <  a  <  1  and  0  <  c  «  1  . 


8.4  A  Partial  Strategy  for  Equilibrated  Reduced  Matrices 

Clearly,  if  A  =  and  each  of  its  reduced  matrices  Is 

equilibrated  (i.e.  max  for  each  1  ),  then  the  partial 

i,j 

strategy,  whereby  we  choose  the  2x2  principal  submatrix  whose  off- 

(k) 

diagonal  element  is  the  maximal  element  in  the  first  column  of  A  , 

is  stable  since,  by  the  equilibration,  such  a  maximal  element  In  the 

(k)  (k) 

first  column  of  A  is  also  a  maximal  element  of  A 

But  we  would  have  to  equilibrate  A  at  the  start,  and  then  equi- 

(k) 

llbrate  each  reduced  matrix  A 

Let  us  now  consider  such  a  partial  strategy  when  A  Is  equilibrated, 
but  we  do  not  equilibrate  the  reduced  matrices. 

8.5  A  Partial  Strategy  for  Unequilibrated  Reduced  MatrlC'is 

Let  A  be  an  n^n  symmetric  equilibrated  non-singular  matrix  with 
max  |Aj^j  1  pQ  for  each  1  .  We  shall  now  consider  the  partial  strategy 

defined  In  §§8.3-8, 4,  but  we  shall  not  equilibrate  the  reduced  matrices 

A(k) 


for  k  <  n  . 


Let  be  the  reduced  matrix  of  order  k  ;  let  -  A  . 

Let  -  max  .  (We  ahall  not  actually  calculate 

l.j 

Let  max  I  "  ®*y*  We  shall  assume  we  have 

i  ^  ^ 

Interchanged  rows  and  columns  so  that  '  Then  let 

^  (k)  lA(k)|  „  A (k)  ,  ,  (k)  .  ,  , (n)  (n) 

X  •  mux  jAji^  I  .  So  X  -■  ^0  ’  '^hlle  X  “  ^^0  “  ^0  ' 

(k)  (k) 

We  shall  use  a  I'^l  pivot  iff  X  ,  Let  m  be  any 

multiplier  (§3.1). 

If  then  we  use  as  a  l^l  pivot,  and 

|mi  ^  <  l/o  ,  while  max  (a^^  +  X^*^Va  ^ 


i.J 


(1  +  1/a)  . 


If  *■  ®  X^^^  then  we  interchange  so  that  "  X 


Let 

v<'^>  = 

|.(k)  .(k) 

lAfi  A22 

-  A<‘'>'| 
A21  1 

Then 

lm|  1 

,(k)  ^^(k) 

+ 

max 

i.J 

iA{;-^>i 

1  -  (k) 

1  IWq  - 

2X^'‘>^  (p<' 

.  So  >  X(’'>^  -  >  (1  -  a^) 


and 


.(k) 


As  In  §5.7  we  would  choose  a  -  “  (1  +  /rT)/fl  .  Then 

max  I  <  (2.57)"  ^  Uq  for  each  A^'*’^  . 

l.j  ^ 


(k) 

But  we  cannot  obtain  a  bound  on  A  as  In  Chapters  11-12  since 

(k)  (k) 

we  cannot  express  bounds  on  V  In  terms  of  . 


^(k) 
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We  leave  the  significance  of  this  method  (§8.5)  In  relation  to 


the  method  In  §§8.1-8. 4  as  an  open  question. 
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Chapter  9  :  Operation  Count 


9 . 1  Solution  by  Diagonal  Pivoting 

Wc  now  consider  the  amount  of  work  required  to  solve  A  X  ■■  B 
by  the  (unequlllbrated)  diagonal  pivoting  method  (Chapter  8),  where 
A  Is  (in  nxfi  nun-singular  symmetric  matrix  and  B  is  an  n^k  matrix, 
l.c.  there  are  k  right  hand  aides. 

In  matrix  notation,  we  perform  the  following  steps: 


(1)  A  =  M  D  . 

(2)  C  =  B  . 

(3)  Y  =  d“^  C  . 

(4)  X  =  Y  , 

— t  *■!  t 

Here  M  means  (M  ) 

Let  p  be  the  number  of  l^l  pivots  used.  So  q  =  -I-  (n  -  p) 
pivots  of  order  2  are  used.  Let  be  the  reduced  matrix  of  order  1  . 

Definition: 


Tl  If  A^^^  uses  a  1*1  pivot 


pivot  [1] 


2  if  A 


(1) 


uses  a  2x2  pivot 


0  if  uses  a  2x2  pivot 


We  shall  use  the  term  Mults  (Adds)  to  mean  the  number  of  multi¬ 
plications  (additions).  We  shall  count  a  comparison  as  an  addition. 

We  shall  use  to  denote  summation  over  those  Indices  1  , 

1  £  1  ;<  n  ,  such  that  pivotli]  "=  j  . 


Summary  of  the  Work  Required 
StepH  (1),  (2),  (3),  (4)  In  59.1  req'ilre: 


9.2 
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Mults  -  ^  +  (k  +  ■|)n*  +  (^  -  2k.)n  +  (k  +  3/2)p  +  31 

1  I  n’  +  (k  +  |■)n*  +  <3  -  2k)n  +  ^(2k  - 

Adds  -  i  +  (k  +  5/8)n’  +  (f  -  2k)n  +  (V  -  3/8)p  +  ^  1  (I  -  1) 


^  3  n’  +  (k  +  3)11^  -  (k  -  I")!!  ~  3  • 

For  k  1  (l.e.  one  right  hand  side),  let  us  compare  the  work 

required  for  the  diagonal  pivoting  method  (applicable  to  all  non-singular 
symmetric  matrices)  with  that  required  for  Cholesky’s  method  (applicable 
only  to  positive-definite  matrices). 


Cholesky 

Diagonal  Plv 

'oting 

Exact 

Lower  Bound 

Upper  Bound 

Mults 

1  3^3  2^1 

T-n  +  —  n  +-rn 
0^3 

i  n^  +  1  n^  +  3  n 

i  +  i  „2  1  n  + 

30  +30  -  3  n  +  32 

Adds 

i 

1  3  .  2  7 

T-n  +  n  “TO 

0  0 

i  n’  +  n'  -  f  n 

4  +  X  -  “T  n 

3  2  0 

1 

Root 

Recipro¬ 

cals 

! 

1 

n 

j 

0 

0 

Since  the  time  required  for  a  computer  to  perform  a  multiplication 
is  much  longer  than  for  an  addition,  the  requirement  of  ~  ^  n^  multi¬ 
plications  and  between  ~  3  n*  and  ~  n^  additions  is  very  satisfac¬ 
tory  for  aymmetri,-  indefinite  matrices. 


IsJit-* 
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Mults 

Mults 


If  A  la  positive  definite,  then  Cholesky’s  method  la  preferable 
to  the  diagonal  pivoting  method.  See  aleo  Appendix  A  (§A.l). 


9.3  Forming  (1)  A  -  M  D 

First  we  decompoae  A  into  A  ■  M  D  . 

Let  us  consider  the  reduced  matrix  A^^^  of  order  1  . 
The  following  chart  shows  our  course  of  action. 


find  :  Adda  -  1  -  1 


find  ;  Adda  -  (1-1) 


■^0 


Mults,  Adds  ■■  1 


yes 


no 


»  1-1: 


use  l^l 

Interchange 

I 

calculate 

multipliers 


,  Adds,  ..  .(l-i) 

’  }  :torm  A 

(1-1) 


"use  2>'2 


Interchange 


1  1  ^  j  ^  rMulta  -  2, 

calculate  det:l.,.  , 

Adds  *  1 

calculate  r  Mults  ■  6(1-2), 
multipliers  ’  Adda  ■  2(1-2) 

1 


1  ;=  i-1 


,  .(i"2)  /  Mults,  Adda 

form  A  :{ 


1  :=  1-2 


Thus  for  A 


Mules  ■  1(1+1)  and  Adds  "  If  plvot[l]  -  1  , 

Mults  "  i*  31  -  7  and  Adds  -  (1-1)  (-j  i+1)  If  pivot  (ij  ■  2 

Recall  that  there  are  p  1x1  pivots  and  q  “  -I-  (n-p)  2x2  pivots. 
For  step  (1): 

Malta  -  j  1  (l+l)  +  (1^  +31-7) 

-  j  1  (1+1)  +  {-  1  (1+1)  +  ~  (1-1)1  +31-7} 

•  l~i  (1+1)  +  (31  -  7)  , 

(9.3.1)  -  ^  n^  +  ^  n’  -  ~  n  +  |  p  +  3  1 

>|n'+|n'+3n 

Adds  -  (1-1)  (I  1+1) 

“  i  i)  +  <i  i'  -  I  + 

{|  1^  +  I  1  +  I  (l-D*  +  \  (1-1)  -  f) 

-  I  +  7  -  I  <«-p>  +  i  i  (^-2) 

i»l 

(9.3.2)  -  in’+|n^-|-n+|p+^  1  (1-2) 

..1  3.5  2  1 

-4  "  +8  "  ■  4  "  • 

Now  let  us  bound  3  i  and  ^  1  (1-2)  from  above. 

1-  1  ?  (p+2J)  -  T-  (n^-p^)  +  ^  (n-p)  and 
3-1 
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t  (i-2>  <  ^  1  (1-2)  -  p{n*-(p+l)n  +  3P*+'|‘P“f^ 

l"n-p+l 

So  we  liave  Che  following  upper  bounds  for  v.1)  : 

"  S  S  3  3 

(9.3.3)  Mules  l^n^+l-n^-^n-^p^-^p  . 

(9.3.4)  Adds  1 n’  +  pH  p  +  |)  n*  -  I  (p*4p+l)n  +  I  ^  p 


For  (2)  C  -  M  H  : 

Multe,  Adds  -  -  n  +  Y  P)  k  . 

(Note  that  1  “  pivot  [1]  »  2.) 

For  (3)  Y  =»  d"^  C  : 

Mults  •«  3n  -  2p  ,  Adds  ■  n  -  p  . 

For  (4)  X  -  Y  : 

Mults,  Adds  «=  (|  -  n  +  ^  p)  k  . 


9.5  Total  Work  Required 

Thus  steps  (1),  (2),  (3),  (4)  together  require: 

Mules  -  ^  n’  +  (k  +  y)  n*  +  (^  -  2k)  n  +  (k  +  y)  p  +  3  1  , 

Mults  £  nH  (k  +  |)  n^  +  (^  -  2k)  n  +  (k  -  |)  p  -  I  p*  , 

Mules  ^  "2^  n^  +  (-j  2k)  n  . 

Adds  -  i  nH  (k  +  |)  n'  +  (|  -  2k)  n  +  (k  -  |)  p  +  i  i  (1-2)  , 

Adds  <  T  n’  +  (k  +  I)  +  (|  -  2k)  n  +  ~  p'  +  (|  "  f)  p'  + 


Adds  >  nH  (k  +  |)  n^  +  (-  -  2k)  n  . 
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9 . 6  Upper  Bound  on  Multe 

We  would  like  a  minimal  upper  bound  on  Multe  that  is  independent 
of  p  ,  where  0  ^  p  ^  n  . 

Let  f(x)  -  vk  -  ^)x  -  J  x==  .  By  elementary  calculus, 

<  f  (2h^)  -  ^^2  -  7)^  .  Thus  Multa  <  |  n*  +  (k  +  |)  n^  + 


(|  -  2k)  n  +  i  n'  +  I  “  I  ^  fl  • 


Upper  Bound  on  Adda 

Now  we  would  also  like  a  minimal  upper  bound  on  Adds  that  Is  Inde¬ 


pendent  of  p 


^  x’  t  -  f)  X*  +  (i  n  -  I  n  +  k  -  ^  )  X  . 


Let  8(x)  -  -j^  X'  t  Ig  -  X  -r  tj  n  -  ^  ^ 

Since  g'<x)  >  0  on  10. nj  we  have 

.  8(p)  £  n  ~  f  ~  12^  "  ‘ 

Adda  <  ^  n^  +  (k  +  |)  n^  -  (k  -  I-)  n 

.in’+-|n^-|n  (where  k  -  1)  . 
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10.1  Introduction 

Let  UB  attempt  to  solve  A  x  ■  b  by  the  diagonal  pivoting 
method.  If  Is  the  solution  we  obtain  from  the  computer,  we  may 

consider  x^  as  the  exact  aolutlon  of  the  system  (A  +  E)  z  -  b  . 

As  In  §2.5,  we  are  Interested  In  showing  that  the  elements  of  E  are 
small  In  comparison  to  the  corresponding  elements  of  A  (Wilkinson, 
1960,  1961,  1963,  1965).  Such  an  analysis  of  E  is  called  a  back¬ 
ward  error  analysis. 

10.2  The  Occurrences  of  Error 

Suppose  that  we  could  perform  the  diagonal  pivoting  method  In 
exact  arithmetic.  In  order  to  solve  A  x  -  b  ,  we  perform  the  fol¬ 
lowing  steps  (see  §9.1): 

(1)  A  -  M  D  (the  decomposition) 

(2)  c“M^b  (the  new  right  hand  side) 

(3)  y"D^c  (solve  the  1><1  and  2x2  systems) 

(4)  X  “  M  ^  y  (recover  the  solution) 

However,  in  finite  precision  arithmetic,  we  have  error  at  each 
step.  Instead  of  decomposing  A  into  MOM*"  ,  we  obtain  M  and  D 
such  that  M  D  A  +  F.  Instead  of  calculating  M  ^  b  ,  D  ^  c  , 

M  *■  y  ,  we  actually  obtain  c  »  (M  +  ^  b  ,  y  ■  (D  +  6  t))~^  c  , 

X  “  (M  +  Mj)  ^  y  for  some  perturbations  ,  6  D  ,  respectively. 
Thus  we  actually  perform: 
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(1)  A  +  F  -  M  1) 

(2)  c  -  (M  +  M^)"^  b 

(3)  y-(P  +  6D) 

(4)  X  -  (M  +  y 

10.3  The  Error  Matrix  E 

Thus  we  have  b  “  (M  +  Mj^)  c  ■  (M  +  Mj^)(D  +  6  D)  y 
-  (M  +  M^)(D  +  6  D)(M  +  x  -  {M  D  (D  +  6  D)(M  + 

+  M  [6  0  (M  +  +  D  1)  X  .  But  M  D  -  A  +  F  ,  while 

(A  +  E)  X  ••  b  .  Henc'i 

E  -  F  +  (D  +  6  D)(M  +  +  M  [iS  0(M  +  M2)^  +  D  ]  . 

If  we  can  bound  the  elements  of  F  ,  M  ,  D  ,  •  and  <5  D  , 

then  we  can  bound  the  elements  of  E  .  We  shall  see  that  most  of  the 
error  lies  in  F  ,  In  other  words,  most  of  the  error  occurs  when 
obtaining  the  decomposition  M  1>  of  A  . 

10.4  Notation 

In  the  following  sections,  the  symbol  e ,0  will  stand  for  any 
numbers  with  |g|  ^  2  ,  |nl  ^  (1.06)  2  ,  where  t  lb  the  num¬ 

ber  of  binary  digits  in  the  computer.  Each  occurrence  of  £  or  n 
in  an  equation  should  be  Indexed,  but  we  shall  suppress  these  indices 
for  the  sake  of  clarity. 

and  shall  denote  the  sumnatlon  over  those  indices  k  , 

1  £  k  ^  n  ,  with  pivot  Ik]  ■  1  and  2  respectively. 


sa 


We  shall  write  g(t)  -  0  (2"*^)  if  11m  2^  g{t)  la  finite. 

t-KO 


10.^  Summary  of  the  Error  Analysis 

In  tlic  following  sections  we  shall  show: 

..(k) 


.  ..Oi) 


(lO.b.l)  |k^j1  2*^  (1  +  3.01/al  Uq  ’  +  [1  +  11.02/(l-a)]  Pq 


for  11',  •  F  -  F  . 


t  /t._\ - X 


(10.5.2)  |e.,|  <  |f  I  +  e  ,  where  e  <  2  (l+o)  max  p 

IJ  IJ  ^ 

max  ^l/c^^  1/(1-Cl^)}  . 


Also,  for  a  “  Oq  ,  we  shall  show; 


(10.5.3)  |f  I  2'-  <  5.71  I;  p^*"^  +  -31.6  l'^  for  1  >  J 


ij 


(10.5.4)  IeI  <  2''-  max  (23.54)n^ 

k 

n 

where  I  I  ie  the  one-norm;  IeI  “  max  '[  1^.  ,|  . 

J  1=1  ^ 

In  Chapters  11  and  12,  we  shall  show  for  a  =  a.  : 


(10.5.5)  max  |F^.|  <  (15.8)n  2  ^  /n  f (n)  (3.07)  n 


0.446 


l.j 


iJ 


(10.5.6)  IfI  <  (15.8)n*  2"*^  /n  f  (n)  Pq  (3.07)  n' 


0.446 


(10.5.7)  >eI  <  (23.54)n*  2“*^  /n  f(n)  p^  (3.07)n 


0.446 


For  Gaussian  elimination  with  complete  pivoting  (see  §2.5): 

(10.5.8)  If!  <  (2.01)n^  2  ^  /n  f  (n)  p^  ,  where  L  U  -  A  +  F  . 


For  further  remarks,  see  §10.12  and  §10.16. 


59 


10,6  The  Decomposition  for  the  Reduced  Matrix 

(r  ) 

Let  A  be  the  reduced  matrix  of  order  r  ,  1  5.  r  5^  n  .  Let 

8  "  plvot[rl  (see  §9.1).  Thus  8  ■  1  or  2  . 

(r ) 

Let  A^  be  the  matrix  resulting  from  deleting  the  first  s 

rowB  and  columns  from  A^*^^  . 


Let  A' 


- \ 

c(r)  1 

8 

8 

S 

e 

,  where  P  Is  s^s  . 
8 


Is  (r-s)>«s  .  We  shall  use  the  s^s  pivot. 


Let  be  the  (r“8)’<8  matrix  resulting  from  calculating 

(.(t)  finite  precision.  Let  be  the  (r“s)5<s 


s  s 


error  matrix 


;  0<*->  -  .  Let 


be  the 


(r )  (r )  (r )  c 

reduced  matrix  of  order  r-s  resulting  from  calculating  A^  ”  *^8  ^s 

( r  ~*  s  ) 

in  finite  precision.  Let  G  be  the  {r-8)>^(r-8)  error  matrix  : 

g(r-8)  .  Jr-s)  _  ^(r)  ^  „(r)  ^(r)t  . 

£  9  & 

But  = 

8  S  S  S  S  8 


S  8  S 


- 

M<^> 

p(r) 

e(r)t 

8 

8 

8 

S 

M(r) 

p(r) 

8 

8 

8 

(r)  p(r) 

e(r)t 

where 


s  s  s 


10 . 7  The  Error  Matrix  for  the  Decomposition 

From  §10.6,  we  have  A+F«MDM^  ,  where  F  ■  and  for 

1  >  J  :  =•  Gjj  +  (M  D  0^)^^  ,  where  for  1  ^  j  : 
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k.2 


i  with 


and  for  j  ^  k  : 


'l-k.J-k 


,(n-k-l) 

'l-k-l,J-k-l 


(g(n-k) 

1 

g(n-k-l). 

'■2  ^j-k-1,1 


,0(n-k). 

'‘2  'j-k-1,2 


If  pivot tn-k+1]  «  1 


If  plvot[n-k+l]  “  2  , 


if  pivot [n-k+lj  ■  0 


if  pivot[n-k+l)  “  1 


if  pivot (n-k+1]  “  2 


if  pivot [n-k+1]  «  0 


10.8  The  Error  Matrlcea  ,  6  D, 


(M  +  Mj^)  c  ■  b  .  From  §10.6, 


'  1  i-j,l 


if  pivotin  j+1]  “  1 


Mij  =  J  ^  ^^\-j-l  1  pivot Jn-  J+l]  -  2 

pivot[n-1+l]  =  0 

V 

for  i  >  j  ,  with  M  .  4  ”  ®  pivotln-J+1]  =  2  . 

»  J 

while  Mjj  *  1  for  every  J  ,  and  "  0  for  J  >  i  . 

Thus  J  1  i  ^**l^j+l  j  “  ® 

pivot tn-j+1]  =  2  . 

(D  +  6  D)y  ”  c  •  From  §10.6,  D  Is  a  block  diagonal  matrix  with 
blocks  of  order  1  and  2.  The  blocV:s  are  the  pivots  in  §10.6, 
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s  "  pivot (r]  .  Also,  <5  D  has  the  same  block  structure  as  D  . 

(M  +  X  «■  y  .  mI)  has  the  same  structure  as  . 

10.9  Floating  Point  Error  Analysis 

Since  most  computers  now  perform  calculations  In  floating  point 
arithmetic  rather  than  in  fixed  point  (see  Wilkinson,  1965,  pp.  110- 
188),  we  shall  give  an  error  analysis  only  in  floating  point. 

Then,  from  Wilkinson  (1965),  pp.  114-117,  we  have  z  =  fX.(x  *  y)  = 
(x  *  y)  (1  +  e)  ,  where  |e|  ^  2  ^  ,  t  is  the  number  of  binary  digits 
used  by  the  computer,  and  *  is  any  arithmetic  operation  (+,  -,  . 

Further,  we  shall  asswoe  that  the  computer  can  accumulate  inner- 
products  in  floating  point  arithmetic.  Then  (fil  (x^^  "*■  ^n^ 


II 

(Xi  yi  +  ...  +  Xn  y^^)|  <  I  i!y^  x^  yj  .  where 


i=l 


IYiI  <  I  (1.06)  .  As  in  §10.4,  we  shall  assume  n  Is  any  number 


-2t 


such  chat  |ri|  £  (1.06)  2  and  we  shall  suppress  indices.  From 
§10.4,  we  may  write  |ril  ^  2  *■  0  (2  ^)  . 


10.10  Floating  Point  Analysis  for  F^ 


consider  A^’^^  .  . 

Case  1 ;  s  *  pivot [r]  ■  1  . 


Ut  1  >  j  . 


and  C^’^^  are  (r-l)xl  .  =  A^^^  = 

So  (1  +  t.)  , 


$ 

I 

/ 


J 
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So  E 


(r-1)  _  /*(r) 


c  -  (2  +  e)  . 


Thus  5.  2'*^  {1  +  (2+2  *^)  }• 

'  (c;*'^  (!+£)•  SO  -  e(c{’^>  . 

i(m;^'>j  i  <1 "  ■  '"I'Ll  ^ 

1  1/a  . 

So  |(Qj’'^)jli  2"Va  and  <  (1  +  2'")/a  .  Thus 

(10.10.1)  <  2"^  fi  +  2/a  +  0(2  *^)]  ,  and 

(10.10.2)  ^  ' 

^(r)  2-t  [i/a  +  0(2"*')]  . 


Case  2:  s  »  plvottr)  -  2 


©2*^^  and  are  (r-2)x2  .  = 


and  |det  P 


^*^>1  -  -  pj’'^  id-  ot")  Pq' 


So  A 


-  t[(A2’'^)ij  -  (1  +  e)  - 
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(M^‘^^12  <C2'^^j2  •  *''^°®  510.6,  we  have 


(10.10.3)  IgJ^'^^I  <  2“*^  {1  +  [3  +  0(2"‘^)]  + 


'Ij  '  -  "  ’^(J 

(2  +  2’‘'l  ^ 


Now  (M<">)  -I  -7V-775 - 7^yr^ -  1  (1  +  e)  . 

Thus  [k\1^  .  IA;;>  A<^2^  -  A<^2^  )  + 

4^  e  (2  +  e)  -  A^2^  A^J^  e  (2  +  c)  +  A^J^^  c  -  A^^^  e  . 

From  §10.5, 

l(0^''^lll  1  (4^  +  2~^  (2  +  2’'')  +  2“*^  +  2'‘^]/v 


2"’^  (2  +  2"'")  2"'^  (1-Kx^)  „-t 

<  - ;— r - ; — r -  + - 7 — ^ t -  2  +  2  '”  +  ,  ^ 

,(r)  ,  (r)  ,,  .2x  .  (r)^  -  1-a  l+a 


>^0  -  ‘^l 


(l-aO  M 


,-t  .  1-Ki" 


But  (1  +  a^)/(l  +  a)  <  1  since  0  <  a  <  1  .  Thus  we  have 


(10.10.4)  l(02'^^)ill  i  2~^  {3/(1  -  a)  +  0(2"")} 


Similarly,  we  obtain 


(10.10.5)  ~  a)  +  0(2"'-)}  . 


(r) 
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From  §10.6,  we  recall:  +  0^'^  .  But 

(10.10.6)  1(0^*"^  I  -^0*^^  ^^0*^^  <  1/(1  -  a)  . 

Thus  from  (10. 10. A)  -  (10.10.6),  we  have  for  j  -  1,2  ; 

(10.10.7)  |(M2‘'^>ijl  1  +  2"‘'  [3  +  2"‘']}/(l  -  a)  . 


From  (10.10.3)  and  (10.10.7),  we  conclude; 


(10.10.8)  <  2'''^  Uq’^^  {1  +  5/(1  -  a)  +  0(2'*^)}  . 


How  we  need  to  bound  ^2^^ 

recall:  ■=  -F  0^*'^  P^*^^  . 


From  §10.6,  we 


From  (10. 10. A),  (10.10.5),  and  §10.6,  we  have 

(10.10.9)  ^  ^  -  a)  +  0(2"*')]}  , 

(r )  (t ) 

since  ^  <  a  \x^  ^  . 

From  (10. 10. A),  (10.10.5),  and  (10.10.9),  we  conclude: 

(10.10.10)  |(M^*'^  P^*^^  qJ*^^*')^^!  <  l3  +  0(2'*')}  2  y^*'^  2"*'/(l  -  a)  . 

10.11  Summary  of  Floating  Point  Analysis  for  F 


From  §10,7,  (10.10.8),  and  (10.10.10)  we  have  for  1  ^  .1  : 

(10.11.1)  iG^jl  2"  <  [1  +  2/a  +  0(2"*'))  l'^  + 

tl  +  5/(1  -  a)  +  0(2"*')}  I'l  , 

(10.11.2)  1(M  D  2*'  <  [1/a  +  0(2"*')]  + 

[6/(1  -  a)  +  0(2"*')]  I'l  . 
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From  (10.11.1)  and  (10.11.2),  we  conclude  for  i  i  : 

(10.11.3)  iF^^l  2*^  <  U  +  3/a  +  0(2~‘^)1  + 

[1  +  11/(1  -  a)  +  0(2'*")]  I'l  . 

We  shall  now  aup-ime  tha;.  the  first  0(2  ^)  term  in  (10.11.3) 
is  bounded  by  0.01  (which  is  true  for  t  ^  8  for  all  a  )  and 
that  the  second  0(2  *■)  term  in  (10.11.3)  is  bounded  by  0.02  (which 
is  true  for  t  ^  14  when  a  ■  (1  +  *Tf/8)). 

Under  these  assumptions,  we  have: 

(10.11.4)  If^jI  2*^  <  (1  +  (3.01)/al  l'^  + 

[1  +  (11.02)/(1  -  a)]  J;;  , 

while  for  a  »  ,  we  have: 

(10.11.5)  If^jI  2^  <  5.71  +  3.16  . 

10.12  Comments  on  the  Bound  for  F^ 

From  (10.11.4),  (10.11.5),  and  §10.3,  we  see  that  we  have  a  good 

(Ic) 

bound  on  F  (and  hence  on  E  )  only  if  the  roaxioial  elements  in 

(k) 

the  reduced  matrices  .>  do  not  grow  too  large. 

From  Corollary  3  of  §5.7,  we  have  for  a  «  : 

<  (2.57)”  Pq  .  Thus  (10.11.5)  does  not  automatically  give  a 

good  bound  on  F  if  n  ^  t  .  At  the  present  time,  most  computers 
have  t  V  50  .  Hence  (10.11.5)  cannot  guarantee  a  good  bound  on  F 
for  systems  of  order  ^  50  . 
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According  to  Wllltlnson  (1965,  p.  215),  If  L  U  -  A  +  F  ,  then 

1^1  J  -  ^2.01)  0  2*"  max  for  Gaussian  elimination,  where 

l.j  J  k 

m^x  ^  2”  for  partial  pivoting,  while  ^  £  I'lt  f(n)  Pq 

for  complete  pivoting. 

In  Chapter  11,  we  shall  show  that  for  the  diagonal  pivoting  method 
with  0  <  a  <  1,  we  have: 

max  p^*^^  <  </n  f(n)  p.  c(a)  h(n,a)  . 
k 

In  Chapter  12,  we  shall  show  that  for  a  “  : 

0(0^)  h(n,aQ)  <  3.07  (n-l)^'^^^  for  n  >  2  . 

Then  from  (10.11.5)  we  have  for  a  “  Oq  : 

l^n!  1  2“^  f(n)  Pq  (3.C7)n®'^^^  [5.71p  +  31.6(n-p)/2]  , 
l.j  ^ 

where  p  is  the  number  of  l^l  pivots  used.  So 

“ax  If^J  1  2'^  /iT  ((n)  Pq  (3.07)n°'^^^  (15. 8n  -  10.09p] 

l.j  ^ 

£  15. 8n  2’'^  /n  f  (n)  Pq  (3.07)n°'^*^  , 

0.446  I  I 

which  is  within  a  factor  7.9(3.07)n  ’  of  the  bound  for  max  |F  ,| 

l.j  ^ 

for  Gaussian  elimination  with  complete  pivoting. 

10.13  Floating  Point  Analysis  for  ^ 

The  blocks  of  D  are  of  order  1  or  2  . 

(a)  Suppose  is  a  block  of  order  1.  Then  so  is  ((6  0)^^^)  . 

We  want  to  solve  ^1  ”  *"1  ’  ol*l*lt>  ^^ii  ^  ^^ll^^i  *  ^1 

due  to  the  finite  precision;  plvotIn-1+1]  •>  1  and  . 
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Thu8  (1  +  e)  and  (6  D)  11  ’  -  “ll  ^  • 

So  1(5  D)^J  <  2’*^  |d^J  .  But  IDjJ  -  ^(n-i+1)  <  . 


Thus  I  (5  D)^J  <  2'*^  . 


(b)  Suppose 


11 


1+1,1 


*^1+1,1  “l+l,i+l 


p(n-l+l)  ^  block  o£  order 


2.  Then  pivot (n- 1+1]  •  2  and  ^  <  a  p^ 

±  Vl,l+1  <1  -  ‘'l+l,!  ^>1  (1  +  O  . 

where  n  -  2“*^  0(2"*^)  .  Since  je  |<  2~^  ,  we  have 

|^<n-k+l)  _  <  2'*'  tl  +  0(2'*')]  (1  +  a^)  . 

.(n-k+1)  j 


Now 


"ll  (1  *  ">  -  'l  Vl.l  “  *  "OO  * 


•"<*  Vi  -  I'j  Vi.ltl"*''*  -  ‘i*i  ■>1*1.1  <■*’'”  ■ 

since  |el  <  2“*^  ,  ln|  1  2"*^  0(2“*^)  ,  and  P^"  <  a  p^" 

after  much  manipulation,  we  have: 


1(5  U)iJ 


1(5  D) 


1+1,1 


°\,1+1 

1(5  D) 


1+1,1+!' 


<  2-'  [1+0(2-^)] 


10.14  Floating  Point  for  and 

'  c 

From  §10.2,  (M  +  M^)c  -  b  and  (M  +  M^)  x  -  y  . 
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From  Wilkinson  (1965,  pp.  247-249),  we  have: 

0(2'")!  Nil  I  *  ^“l\j  “  ° 

for  i  <  j  , 

.  iu-.2)jj|  <1  2"*^[1  +  0(2"‘')]  (n  -  2  +  1  -  j)  \h^^\ 
for  1  ■>  J  . 

From  §10.10,  »  1  and  for  i  >  j  : 

Nij  I  1  ®ax  (l/a,  l/(l-a)}  [1  +  0(2~*')]  . 

So  for  i  >  j  : 

(10.14.1)  |(Mp^jl,  |(M^)^j|  <  2"*^  max{l/a,l/(l-a)}[l+0(2"‘')]3(n-2+i-J)/2 
and  |(M^)^J  ,  <  2*''  H  +  0(2“S]  . 


10.15  Floating  Point  Bound  for  _E 


From  §§10.3,  10.11,  10.13,  and  10.14,  wc  could  exprcsr  In 

terms  of  F^^  and  the  elements  of  M,  D,  M2,  and  6  D  .  But  this 
expression  is  very  complicated. 


Let  us  define  IaI 


n 

I  |a,.1 

j  1=1  J 


where  A  Is  n^n 


(This  is 


usually  called  the  one-norm.) 


Clearly,  Ni j I  i  • 


From  §10.3,  we  now  have 
(10.15.1)  |e^j|  <  |f^^1  +  E  ,  where 

e  -  Id  +  6dI  I (m+m^)!  +  ImI  [I6dI  I(M4M2)I+  IdI  Im^I  ]  . 
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lul  <  (1  +  u)  max  •  Hrom  010.13, 

k 

I6dI  2  ^  (1  +  a)  max  II  +  0(2  *■)  ]  . 

k 

Fioiii  £10.10,  ImI  ,  £  1  +  (n-1)  max  {l/a,  l/(l-a)}  (1  +  0(2  *■)]  . 

Let  US  assume  0(2  *")  (1  +  max  (l/ot,  l(l-0()}  3(n  +  2)^>«  0(2  , 

i.e.  let  us  assume  n  2  <<  1  .  Then 

IMjI,  Im'I  <  2''*^  (1  +  0(2“'^)]. 

Now  we  can  bound  e  .  Let  6  •  1  +  (n-1)  max  (l/a,  1/ (1-a))  . 

(10.15.2)  G  <  2”*^  (1  +  a)  max  6(2  +  6)  (1  +  0(2"*^)]  . 

k 

For  a  •  ,  we  have  6  <  2.781n  ^  ao 

(10.15.3)  e  <  2"*"  max  (7.734)n*  U  +  0(2**')]  . 


i 

! 


3 

i 

I 

I 
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10.16  Comments  on  the  Bound  for  E 


From  §10.12  and  §10.15,  we  need  to  bound  the  In  order  to 

have  a  good  bound  or  E  . 

(k)  r~ 

In  Chapter  11  we  shall  show  that  max  Pq  1  /n  f (n)  p^  c(a)  h(n,a) 

k 

for  the  diagonal  pivoting  method  for  0  <  a  1  .  In  Chapter  12  we 

0  446 

shall  show  that  cla^)  11(0,01^)  <  3.07(n-l)  ’  .  Then  we  would  have 


Itl  ■  (23.54)n^  2  ‘  f(n)  p.  {3.C7)(n-l) 


0 . 4  46 


For  Gaussian  elimination,  In  order  to  solve  A  x  “  b  ,  In  finite 
precision  we  actually  perform: 

(1)  L  U  -  A  +  F 

(2)  c  -  (L  +  5L)“^  b 

(3)  X  =  (U  +  5U)~^  c. 

From  Wilkinson  (1965,  p.  215,  pp.  248-252),  we  have  (cf .  §2.5): 

I  —  2.01n  2  max  p'  and 

l.j  ^  ^ 

IeI  <  2.01n*  2"*^  max  [1  +  0(2”'))  . 

k 

(k) 

For  complete  pivoting,  max  p^  £  rn  f (n)  Pq  (Wilkinson,  1961). 

k 

Thus  our  bound  on  IeI  for  diagonal  pivoting  differs  by  a  factor 

36(n-l)*^’^^^  from  the  bound  on  IeI  for  Gaussian  elimination  with 
complete  pivoting. 
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Chapter  11  :  Art  A  Posteriori  Bound  on  Element 
Crowth  fot  0  <  Cl  i  1 

11.1  Introduction 

We  saw  In  Chapter  10  that  the  bouna  on  the  elementa  of  the  error 
matrix  depends  on  the  maximum  of  the  elements  of  the  reduced  matrices. 

^  tl ) 

Let  A  =  A'  be  the  original  aynuietric  matrix  of  order  n  v/ich 
(k) 

det  A  1*  0  .  Let  A  be  the  reduced  nsacrix  of  order  k  . 

We  saw  in  §8.4  that  max  1  <(2.57)”  ^  max  jA.  ,|  .  Cut  we 

i,j  '  i,;)  ^ 

shall  show  in  5S11.2-11.7  that  we  can  get  a  much  better  bound  by  the 
use  of  Wilkinson's  techniques  for  Gaussian  elimination  with  complete 
pivoting  (Wilkinson,  1961;  pp.  281-285). 

His  proof  depended  on  the  fact  that  the  pivots  in  Gaussian  eli¬ 
mination  with  complete  pivoting  were  maximal  elements  in  the  reduced 
matrices.  Our  pivots  are  not  necessarily  maximal  elements,  but  they 
are  closely  related  to  the  maximal  elements  in  the  reduced  matrices. 

We  shall  assume  that  we  use  strategy  (§5.7)  for  any  a  , 

0  <  a  <  1  ,  and  that  for  all  (see 

Lemma  4,  §5,5).  In  particular,  this  lower  bound  holds  for  the  strate¬ 
gies  In  56.1,  §6.4,  .ci.d  §8.1. 

11.2  The  Pivots 

(k) 

Let  A  be  the  reduced  matrix  of  order  k  . 
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,  (k)  I ,  (k.)  I  (k)  I  (k)  I 

Let  vC  “  MX  |a'  I  and  y'  -  max  |a'  |  . 

a  ^  j  ij  1  11 

lat  ufl  assume  that,  whenever  we  shall  use  a  2>‘2,  Interchanges  have 
ensured  that  u  >=  -  A^^  I  i 

(In  particular,  this  assumption  holds  Cor  the  strategies  In 
556.1-6.2,  556.4-6.5,  and  §58.1-8.2.) 

From  strategy  (§5.7)  and  §9.1,  we  recall: 


plvct(k] 


Let  us  now  define 
(k) 


r.i 

) 

i£ 

pf> 

< 

a 

(k) 

i.e. 

if 

uses  a 

l^^l  pivot 

'I  2 

if 

< 

a 

i.e. 

if 

A<‘‘> 

uses  a 

2x2  pivot 

lo 

If 

(k+U 

1 

< 

a 

kr“  • 

i.e. 

if 

^.(k-tl) 

uses  a 

2x2  pivot. 

If  pivot [k]  ■  1 


if  pivot (kj  “  2 
pivot  Ik]  -  0 

The  will  be  called  pivots. 

From  the  decomposition  A  M  D  we  see  that  |det  a|  ■  p^. 
and  jdet  A^^^  1  »  pj^...Pj^  ,  since  det  M  »  1  . 

11.3  Hadamard's  Ineq\;ality  -- 

B>  Hadamard's  Inequality  (Gantiaacher ,  pp .  253-254), 


(11. 3.1)  jdet  a|  <  {  n  I  <  (n  -  {/K  . 

1-1  J-1  y  u 

since  *  'wg  "  |  .  Also 

C ».] 

(11.3.2)  IdfetA*-*^^!  <  vs:  . 


.  (n)vn 


.  Also 
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11.4  Bounding  det  A 

(1)  If  plvotlkl  -  1  :  >a  •  So  Uq*^^  i  p,^/a  . 

(11.4.1)  pj  ...  p^  -  |det  <  (A  <  (p^  ^/a)*"  . 

(2)  If  pivotlkj  ■  2  :  ^  asoume 

(k)  (k)^  (k)  ,  ,,  ,(k)  .  (k)  ,  ,,  2»  (k)^ 

_  ^^0  ”  ^1  ^  1  (1  "  o  )  Wq 

Thus  <  ^^*"^(1  -  a^)  -  p^/(l  -  a^)  . 

(11.4.2)  ...  pj^  =  Idet  <  (p^  A/a  -  a^)  )*"  . 

11.5  Fundamental  Inequality 


From  §11.4  (and  later  in  §11.5),  we  are  led  to  the  following: 


Definition: 


(11.5.1)  6, 


1/(1-Qt*) 


if  pivot tk)  •  1 
if  plvotlkl  ■  2  . 


1  fk+1  il+l/k  ,  f. 

k  plvotlkl  -  0 


From  (11.4.1)  and  (11.4.2)  we  have: 

(11.5.2)  pj^  . ..  Pj^  <  (f^  Pj^)*^  If  plvotlkl  j*  0,  1  <  k  <  n  . 

We  would  like  to  have  a  similar  equation  for  plvotlkl  -  0  for 
ouv  analysis  in  §11.6. 

If  plvotlkl  -  0  ,  then  pivot Ik+ll  -  2  .  By  (11,5.2), 


‘’k  Pk+1  - 


^+1  Pk+1^  •  ^k+1  -  '  a') 


Pk  ■  Pk+1  •  ^  ^  t’y  (11.5.1), 


“k  k 
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■>1  Pk  i  K"*!)  "1:  ■  "k’"  ■  ">"• 


(11.5.3)  ...  P|^  1  (A  k  ,  1  £  k  £  n  . 


11.6  Bounding  Pivot  Growth 

Define  •=  log  .  From  (11.5.3), 
k-1 

(11.6.1)  I  q  <  (k-1)  q.  +  ^  k  log  (k  6.  )  . 

l«il  ^  K.  Z  IC 


Since  |det  •••  ,  we  have: 


(11.6.2)  log  |det  A 


(n). 


n 

I 

1-1 


Dividing  (11.0.1)  by  k(k-l)  for  2  £  k  <  n-1  and  (11.6.2)  by  n-1 
and  adding  we  have ; 

qi  +  *  ^2^^  •+••••  q^_2/(n-2)  +  q^_j^/(n-l)  +  q^/(n-l) 

-  'ih  ^  2 


+  q,/2  +  q.j/3  +  ...  +  q^_j^/(n-l)  , 


n-1  .  - 

after  observing  that  7  —7 — 77-  ■  7 — ; - r  . 

r-k 


n-1 


(11.6.3)  +  q^/(n-l)  £  i  log  n  (k  Bj^)^ 

k»2 

From  (11.5.3),  |det  A^"^|  <  (/rTl”  p  )"  . 

n  n 

(11.6.4)  f(r)  -{  n  , 

k-2 

(11.6.5)  h(r,a)  -{  R  g  l/A-l)}l/2 

k=»2  ^ 


+  log  |det  A 


(n) 


We  define 
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From  (11.6.3)  we  now  have 

<lt  +  q  /  (n-1)  ^  log  f(n-l)  +  log  h(n'l,a)  +  . .  "  , log  (n  0  )  +  n  q  /(n-1)  . 

With  simple  manipulation  we  have 
qj^  “  1  log  f(n-l)  +  2 (n-1)  h(n-l,a) 

®n  I 

-  log  f(n)  +  log  h(n,a)  +  \  log  (n  0  )  . 

i  n 

Thus  we  conclude 

(11.6.6)  p,/p  <  £(n)  h(n,a)  . 

in  —  n 

Similarly,  we  have  for  1  <  k  £  n  ; 

(11.6.7)  Pj^/P^  1  »^n  -  k  +  1  f (n  -  k  +  1)  h(n  -  k  +  l,a) 

£  f(n)  h(n,a)  . 

(11.6.6)  and  (11.6.7)  hold  for  all  a  ,  0  <  a  <  1  ,  under  strategy 
(§5.7)  provided  that  for  all  (§5.5). 

We  now  have  a  bound  on  pivot  growth.  But  we  are  Interested  In 

bounding  element  growth  (§11.1).  If  A  finishes  with  a  2x2  pivot  we 

(2) 

are  Interested  In  bounding  Mq  ,  while  If  A  finishes  with  a  ixl  pivot, 
we  must  bound  . 


11.7  Bout.din,.  Element  Growth 

We  shall  now  express  p, ,  p„,  p  ,  and  T  p  In  terms  of 

1  i  n  n  n-i  n 


/I)  /2) 

0  *  ^0  ' 


0 
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Ut  A  -  ,  Uq  -  .  V  - 


Since  pivot  [n]  0  ,  p  ■  \ 


Uj^  If  pivot [n]  -  1 


/j  if  pivot [n]  ■  2 


,  and 


8n 


1/a*  If  plvotto]  ■  1 

1/(1  -  a*)  if  pivottnl  -  2 
But  £  Mq  always.  If  pivot [n]  ■>  2,  then  <  a  Up  ,  and 
(by  §5.5)  /j  <  /Uq  +  yj  1  ''I  +  a*  Mq* 


Pn- 


If  pivotln]  -  1 


/I  +  a*  Pp  if  plvotto]  -  2 


L 


l/a 


Now  1/(1  -  ot*) 


if  pivot tn-1]  ■  1 
if  pivot [n-l]  •  2 

(n  if  plvotIn-1]  -  0 

n-l 


Let  U6  define: 

(11.7.1)  m(a)  -  max  ,  /(l+a*)/ (1-a*)  }. 
Thus  we  have: 

(11.7.2)  P„  i  m(a)  Pp  ,  and 

n  ^  V 

m(a)  Prt 


(U.7.3)  .^P„i 


if  pivot [n]  ■  1 


^“"^Vn/  (n-l)  pf+STTCr^Poif  plvot(n]  -  2 


From  (11.6.4)  and  (11.6.5),  we  have 
(11.7.4)  £(n-l)  h(n-l,a)  -  f  (n)  h(n-l,a)  . 
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We  have  two  caaes  (§11.6)  as  to  whether: 

(a)  the  last  pivot  is  1^1,  i.e.  plvot[l]  "It 

(b)  Che  last  pivot  Is  2^2,  i.e.  pivoc[2]  ■  2  . 

Case  (a)  :  Let  pivottl]  ■  1  .  Then  pj^  ■  . 

From  (11.6.6)  and  (11.7.2)  we  have 

(11.7.5)  yj^^  -  1  •'JT  f  (n)  Wq  «(a)  h(n,a)  . 

Case  (b)  :  Let  pivot [2]  •  2  .  Then  p^  “  and 

y^^^  <  a  y^^^  .  Since  we  have  assumed  ^  y^^^  -  y^^^  , 

>  (1  -  a*)  yj^^*  .  Hence 

(11.7.6)  1  P2/''l  -  • 

From  (11.6.7)  with  k  ■  2  and  from  (11.7.6)  we  have 

(11.7.7)  y^^^  <  f(n-l)  h(n-l,a)  >T~T  p  /A  -  . 

u  ~  n**i.  n 

lk>w  we  have  two  cases  as  to  whether  Che  first  pivot  Is  (1)  l^l  , 
or  (11)  2x2. 

(I)  Let  plvot[n]  ■  1  •  From  (11.7.3)  and  (11.7.), 

£  /n-l  f(n-l)  h(n-l,a)  m(a)  y^/Zl  -  a*  . 

But  f(n-l)  <  f(n)  and  h(n-l,a)  <  h(n,a)  .  Thus 

(11.7.8)  yj^^  <  f(n)  y^  m(a)  h(n.a)//l  -  . 

(II)  Let  pivot[n]  ■  2  .  From  (11.7.3),  (11.7.4),  and  (11.7.7), 
£  /n  f (n)  Pq  h(r,a)  A  +  a^/(l  -  a^)  . 
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But  /I  +  a^/ (1  -  a*)  £  m(a)//l  -  a'^  from  (11.7,1).  Hence 
(11.7.9)  f(n)  Wq  »(a)  h(n,a)//r^^  . 


Let  us  now  define: 


(11.7.10)  c(a)  -  m(a)  x. 


If  pivot [1]  -  1 


l//r^^  If  plvot(2]  -  2 

m 

where  m(a)  •  max  ll/a,  J (14a^ )/ (1-a* )  }  . 


Then  we  conclude  for  0  <  a  <  1  : 


(11.7,11)- 


If  plvotll]  ■  1  :  y 


(1) 


I  /o\  i  £  ftn)  Oft  c(a)  h(n,a)  . 

1  If  plvotl2l  -  2  ;  Uq  M  ° 

(k) 

Similarly,  we  have  for  each  reduced  matrix  A  , 
(11,7.12)  £  /n  -  k+j  f(n  -  k+j)  jJq  c(k,a)  h(n  -  k+j ,  a)  , 

1  if  J  -  1 


where  J  «  pivot {kj  ,  c(k,a)  “  m(a)  x 


1//1  -  a*  if  j  -  2 


and  0  <  a  <  I  . 


,<k) 


Hence,  for  all  A  ,  we  have  for  0  <  a  <■  1  : 

(11.7.13)  <_  i/n  f(n)  Og  c(a)  h(n,a)  . 


11 . 8  Conments  on  the  Bound 

The  bound  in  (11.7.13)  holds  under  strategy  ,  0  <  a  <  1  ,  and 
(k)  (kl^  (k)*  fkl 

under  the  assumption  v'  i  ~  ^1  ‘  parti¬ 

cular,  (11.7.13)  holds  for  the  unequilibrated  diagonal  pivoting  stra¬ 
tegy  (§§8.1-8. 2),  the  complete  strategy  (§§6.1-6. 3),  and  the  partial 
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tiqullibracud  Mtrntegy  (9S6.4-^.3).  (For  the  partial  equilibrated  stra¬ 
tegy,  we  asaunc  that,  given  a  reduced  matrix  with  ■  mx  |aJV  |  , 

l.j 

(k) 

we  equilibrated  A  so  that  the  maximal  element  in  absolute  value  In 

.  .  (k)  X 

each  row  is  .) 

WllWlnson  (1961)  obtains  1 (n)  aa  the  bound  on  the  elements 

in  the  reduced  matrices  for  solving  A  x  -  b  ,  p  -  max  1a  ]  ,  by 

l,j 

Gaussian  elimination  with  co<s^plete  pivoting.  We  have  the  extra  factor 
c(a)  h(n,a)  since  our  pivots  are  not  necessarily  maximal  elements  of 
the  reduced  matrices. 

We  call  the  bounds  in  (11.7.11)  -  (11.7,13)  a  posteriori,  since 
we  cannot  calculate  the  terms  of  h(n,a)  until  we  know  the  posi¬ 

tion  of  the  blocks  of  order  1  and  2  in  D  for  the  decomposition 
A  ■  M  D  .  In  Chapter  12,  we  shall  give  a  bound  on  c(a)  h(n,a)  , 
independent  of  the  structure  of  D  ,  for  the  value  a  ■  - 

(1  +  /r7)/8  (§5.7). 

11.9  Smaller  Bound  on  Pivot  Growth 


If  H  is  an  n^n  positive-definite  Hermitian  matrix  and 
is  the  minimum  eigenvalue  of  A  ,  then  (Shlsha ,  p.  173): 

(ll.S.l)  det  H  <  n  H  -  (  I  |H,,r)  . 

i-1  niJ>i>l 

If  A  is  an  nxn  non-singular  real  symmetric  matrix  and  X  is 
the  miniaium  of  the  absolute  values  of  the  eigenvalues  of  A  ,  then, 
setting  H  “  A*'  A  in  the  above,  we  have 
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(11.8.2)  |det  aT  i  '1  <  I  I  a  I 

1-1  j-1  J  ny>lll  ^ 


where  Aj  is  the  column  of  A  . 


Then  |det  a|^  £  (n  u^)  e (A)  ,  where 


(11.8.3)  g(A)  -  1  -  X 


2(n-2) 


L 


I  lA 

n^>l>ll 


5*1 1 3 


/(n  Wq)" 


Using  an  analysis  similar  to  that  in  §§11.2—11.7,  we  obtain: 


(11.8.4)  ' 


If  plvotll]  -  1  : 


If  plvotl2]  -  2  ; 


£  /n  f(n)  Mq  c(a)  h(n,a)  t(A)  , 


where  t(A)*  -  e(A)""^  H  e(A^'^^ and 


r-2 


e(A 


(r) 


)  -  1  -  I  ^  ’  where 


(r) 


is  the  reduced  matrix  of  order  r  (A  -  A^**^)  ,  X^  la  the  minimuni  of 
the  absolute  values  of  the  eigenvalues  of  A^  ^  ,  and  Aj  ^  is  the 
column  of  A^*^^  . 

If  det  A  )*  0  ,  then  |X^|  >  0  for  each  A^*^^  ,  If  A^*^^  Is 

not  lladamard  (see  §12.6),  then  e(A^'^^)  <  1  .  If  A^  ^  is  Hadamard, 
then  A^*^^  will  use  a  1^1  pivot  and  will  not  be  Hadamard 

(see  Appendix  A) . 

Thus  t(A)  <  1  if  det  A  >*  0  .  Hence  (11.8.4)  gives  a  lower  bound 


than  does  (11.7.11),  but  (11.8.4)  is  so  complicated  we  are  unable  to  use 
it  to  advantage,  and  we  present  it  merely  for  its  academic  interest. 


similarly,  for  Gauaelan  eltolnatlon  with  complete  pivoting,  we 
can  obtain  £ (n)  t(A)  aa  the  bound  on  the  elements  In  all  the 
reduced  matrices  when  solving  A  x  •  b  ,  V>.  >  max  |a  |  ,  det  A  0 

“  j  4  ij 


if  we  replace  the  In  t(A)  by  jo^| 


of  A 


(r) 


the  singular  value 


of  minimum  modulus. 
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Chapter  12  ;  An  A  Priori  “Vaund  on  Element  Growth 
for  a  -  Oq  -  (1  +  /iT)/8 


12.1  Introduction 

In  {11.7.13)  we  obtained  /a  £(n)  c(a)  h(n,a)  as  a  bound  on 
the  element  growth  In  the  reduced  matrices,  for  any  a  with  0  <  a  <  1  . 

From  (11.5.1),  (11.6.4),  111. 6. 5),  (11.7.1),  and  (11.7.10)  we 
recall : 


£(n)^ 


1/  (f 

1/(1  -  a^) 


1  ,  k+1  ,1+1/1<- 

k  ^ 

k-2 


if  pivotik]  ■  1 
if  pivot(k)  -  2  , 

if  pivot(k]  ■  0 

h(n,a)*  -  n 
k-2 


» 


1  if  pivotlij  -  1 

m(a)  -  max  }  ,  and  c(a)  -  n)(a) 


1^/ /l  -  a'^  if  pivot  [2]  -  2 

The  term  c(a)  h(n,a)  arose  from  the  fact  that  our  pivots  are  not 
necessarily  the  maximal  elements  of  the  reduced  maf'lces,  but  can  be 
expressed  as  multiples  (Involving  a)  of  such  maximal  elements.  The 
bound  in  (11.7.13)  Is  a  posteriori  since  we  cannot  calculate 
c(a)  h(n,a)  until  we  know  the  pivot  selection,  i.e.  until  we  know  the 
position  of  the  blocks  of  order  I  and  2  ■*  a  the  block  diagonal  me  trlx  D 
for  our  decomposition  A  -  M  D  for  a  given  value  of  a  , 


0  <  a  <  1  . 
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We  would  like  an  a  priori  bound  on  the  element  growth,  i.e.  a 
bound  independent  of  the  selection  of  a  1><1  or  a  2^2  pivot  at  each 
stage. 

We  would  also  like  the  a  priori  bound  to  hold  for  all  a  , 

0  <  Gt  <  1  .  But  no  such  bound  exists,  as  we  shall  now  show, 

Let  p  be  the  number  of  I’^l  pivoto  for  the  deccmpi<sltion  (§9.1), 
I.e.  there  are  p  blocks  of  order  1  in  D  above.  If  p  <■  n  (e.g. 
for  a  positive  definite  matrix,  see  .^pendix  A),  then 

h(n,a)  “  n  -v  «  as  o  0  .  If  p  »  0  (see  S&.d),  then 

k-2 

h(n,ci)  >  n  (l//l  -  «  as  a  1  . 

k-2 

Thus  we  shall  give  on  a  priori  bound  for  element  growth  only  for 
the  value  a  -  -  {]  +  i^7)/£  (§5.7). 

12.2  Lower  Bound  on  c(a)  h(n.a)  for  0  <  a  <  1 

We  shall  now  find  a  lower  bound  on  c(a)  h(n,U)  for  all  (X  , 

0  <  a  <  1  ,  in  order  to  show  that  the  upper  bound  for  a  -  that 
we  obtain  in  §12.4  is  a  reasonable  bound,  i.e.  chat  some  other  choice 
of  a  would  not  provide  us  with  a  much  better  upper  bound  on 
c(a)  h(n,a)  . 

Since  min  c(a)  h(n,G()  ^  min  c(a)  ^  min  h(n,a)  ,  we  need 
0<a<l  0<a<l  0<a<l 

only  find  lower  bounds  for  the  latter 

Lemma  1;  min  c(a)  >  2.029  . 

0<a‘'-l 


two  minima. 


1  +  3//5'  >  2.029  . 


q.  e.d. 


Proof :  rain  c(a)  ■  ci//l  -  ^ 

0<a<  I 

We  shall  later  need  the  following: 

n+1  ,  ,  n  . 

I^enuna  2:  [  ^  <  log  <  I  T  for  n  ^  lo  ^  1  . 

k-itt+l  k“a 


Proof  :  Uy  elementary  calculus,  x  ^ 


k=«a 


kH 


—  dx  ■  log  I  ^  n.m  >  1 


q.e.d. 


How  we  cau  find  r.  lower  bound  for  rain  h(n,ci). 

0<a<l 

.  4  u/  V  l08  •'2  .  0.3465 

Lcicma  3:  min  h{n,0i)  >  n  >  n 

0<a<l 

n 

Proof ;  Let  w(n)  «  \  l/(k-l)  .  If  p  »  n  ,  then 

k*2 


h(a,a) 

n/2  1 

where 

t(n)  •  Jl  ' 

Thus 

h(n.cx>  >  b<a) 

0  -  t(n) 

l/(2j-2) 


>  1  . 


rain  b(oi)  •=  ,  and  is  attained  by  a  ■  1/V2  .  From  Lemma  2, 

0<a<i 

w(n)  >  log(n)  .  Eence  min  h(n,a)  >  ^ 

q.e.d. 

From  Lemmas  1  and  3  we  obtain  a  lower  hourtd  for  c(a)  h(n,a) 
for  0  <  a  <  1  . 

0  346*' 

Theorem  1:  min  c(a)  h(n,a)  >  (2.029)n 
0<a<l 
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.0.446 

In  §12.4  we  chall  show  chat  cCa^)  h(n,aQ)  <  3.07(n-l) 

Thus  our  choice  of  «  “  “q  (S5.7)  will  provide  us  with  an  upper  bound 
on  c(Oq)  h(n,0iQ)  which  does  not  differ  much  from  the  minimum  of 

c(a)  h(n,a)  for  Q  <  a  <  1  . 


12.3  Remarks  on  an  Upper  Bound  for  h(n.a) 

We  cannot  evaluate  the  until  we  know  the  pivotal  selection, 

but  If  we  could  bound  all  the  Independent  of  the  pivotal  selection, 

then  we  can  obtain  an  upper  bound. 


Theorem  2:  If  6,^  i  «  for  all  k^r 
is  independent  of  n,  then  for  0  <  a  <  1  : 

h(n,a)/h<r-l,a)  <  Y(n-l/®*  ,  where  y  - 


,  where  c  >  0 


(r-2) 


-  log  /c 


and 


r  >  3 


n 

Proof;  From  Lemaa  2  of  §12.2,  ^  l/(k“l)  *  log  (  ^ 

k-r>3 


1  ,  ,  n-1  . 

"  l/ra-11,  J/'’  '•  r-2  ^ 

Thus  h(n,a)/h(r-l,a}  <{llc  }  “<0 

k«r 


But 


log  v  log  X 


ior 


>  0 


q. e.d. 


This  shows  that  if  we  can  bound  independent  of  the  pivotal 

strategy  for  k  ^  r  and  If  we  can  consider  the  worst  possible  cases 
of  pivot  selection  for  k  ••  2,...,r-l  and  bound  h{r-l,a)  ,  then  we 
have  a  bound  for  h(n,a)  .  (In  order  to  consider  the  cases  for 
2  <  U  <  r-1  ,  we  must  have  r  reasonably  small.) 

In  §12.4  we  shall  do  this  for  a  -  a 


,  and  we  shall  have  r  ■  5  . 
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12.4  ^  A  Priori  Bound  for  a  ■  -  (1  +  >T7)/8 

For  a  -  Oq  ,  max  (1/oIq  ,  /tl  +' 0^'/ (!l '  -  a^))  -  i/a^  .  Thua,  for 

“  "  “o  • 

1  If  plvot[lJ  ••  1 

(12.2.1)  d(aQ)/aQ  ,  where  d(a)  "•<| 

^1/A  -  If  pivot i2]  “  2 

Now  we  shall  show  that  for  a  «»  Oq  ,  is  maximal  for 

1x1  pivots  for  k  ^  5  .  Recall  §9.1  or  §11.2  for  the  definition  of 
plvotjk]  . 


Leiuma 


4:  For  a  *  Oq  and  k  ^  5  ,  if  pivotlk]  »  0  ,  then 


„l/(k-l)  .1/k  ^  .  2.1/(k-l)  +  1/k 

\  \+l  - 


Proof:  From  (11.5.1),  If  pivotlk]  «  0  ,  then  for  0  <  a  <  1  , 


Vi  =  6i,  ‘  r  t 


1  ,  il+l/k 

k  k  ‘  1-a* 


neace 


.l/(k-l)  .1/k 
^k  ^k+i 


2.1/ (k-l)  +  1/k  , 


,  k+1 , .  k 


e,*' eriT  ^  a/ct^y iff  (k+D^'^/k"  <  a^(l/a^  -  1) 


2k 


Iff  (k+1)  log  (k+1)  -  k  log  k  £  log  +  2  k  log  (1/a*  -  1) 
iff  g(k,a)  ^  0  for  k  ^  k^  for  some  Integer  k^  ^  1  , 

where  g(x,a)  •»  2  x  log  (1/a*  -  1)  +  log  a*  -  (x+1)  log  (x+1)  +  x  log  x 


^  g(x,a)  -  2  log  (1/a*  -  1)  -  log  (1  +  1/x)  >  0 

OX 

iff  X  >  l/I(l/a*  -  D*  -  1]  . 

In  order  to  evaluate  this  inequality,  we  now  fix  a  =  .  Let 

G(x)  =  g(x,  Uq)  . 
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Thus  G' (x)  -  ^  ^  ^  ^  ® 

monotone  increasing  function  for  x  >  0.94.  Now  G(4)  ■  -0.476  and 


G(5)  -  3.698.  Thus  gfk.a^)  >0  for  k  ^  5  . 


q.e.d. 


Now  we  can  apply  Theorem  2  of  §12.3  and  Lemma  4  to  obtain  the 
following  two  necessary  lenmaa. 


log  a 


Lemma  5:  If  plvotlS]  ^  2  ,  then  h(n,aQ)/h(4,aQ)  <  3 
<  0.613  (n-l)°'^^^  . 

Proof ;  By  Theorem  2  of  §12.2  and  Lemioa  4,  h(n,aQ)/h(4,a) 

log  a  -log  a  Q 

<3  °  (n-1)  ^  <  0.613  (n-l)*^-^^*"  . 


(n-1) 


-log  a 


q.e.d. 


Lemma  6:  If  pivotlS]  "  2  ,  then  h(n,aQ)/h(3,aQ)  <  0.717  (n-1) 


0.446 


Proof ;  If  pivot(5]  ■  2  ,  then  plvotI6]  2  ,  so  by  Theorem  2  of 

log  a_  -log  a„ 

§12.2  and  Lemma  4,  h(n,a)/h(5,a)  <  4  (n-1)  .  Since 


1  1 


pivot[51  ■  2 


jg,3  g4  j  1/2  <0.717  . 


1  ,  5  ,5/4, 1/3  .  1  ,1/4  ,1/2 


q.e.d. 


2  1/2 

Now  we  need  only  to  bound  {3„  3,  }  '  for  pivot [5]  «  2  and 


1  1 

2  3  1/2 

{62  63  64  ^  for  pivot[5l  ^  2 


Lemma  7:  If  plvotI3J  "  1  ,  then  h(3,ao) 


(I/Oq)  if  pivot [2]  *=  1 


if  pivot  1 2]  =  2 


<  1.96. 


Proof :  Since  plvot[3]  ■  1  ,  S.  ■  ^  and  pivot  I2]  j*  0  . 

“o 


So 


if  plvoCl2]  -  1 


,  and  h(3,a-)  <  1.96. 


^  I-Uq  if  plvot[2]  -  2 


q,  e.d 


Lemma  8:  If  ptvocl3]  •  2  ,  then  hOjO^)  <  2.74. 


Proof:  Since  pivot[3]  -2,6,-  r~r  and  6,  ^  4  ) 

-  3  I-cIq  2  2  I-Oq 


3/2 


q.e.d 


Lemma  9:  If  pivoc[3]  -  0  ,  then 


r 


a/ 3 


h(4.ao)  -  3 


—  If  pivotl2]  ■  1 


if  pivot[2]  *  2 


Proof :  Sines  plvot[3]  -  0,  6^  -  ,  S3  =  ^  ) 


4/3 


and  pivot [2]  j*  0  •  So  ^2  ^  < 


a.. 


i-<x; 


if  pivot 1 2]  =  1 


.  q.e.d. 


_1  if  pivot [2]  ^  2 


Now  we  put  Lemmas  5-9  together  to  get  a  bound  on  cCcIq)  hCn.a^) 


which  ie  independent  of  the  pivotal  selection. 


Theorem  3;  cCa^)  h(n.aQ)  <  3.07  <  3  /JT  for  n  >  2  . 

Proof  5  From  (12.2.1)  •  cCqi^“  * 


d(ao) 


if  ptvotll]  •  1 


1//1  -  If  plvotI2l  -  2 

Now  we  mu8t  conalder  varioua  sltuatlona. 

Case  I:  plvot[5l  -  2  :  Then  plvotP]  0  .  By  Lemma  6. 

.0.4A6 

h(n.aQ)/h(3,aQ)  <  0.717  (n-1) 


(a)  If  pivot [31  ■  1*. 


Then,  by  Lemma  7.  h(3,0iQ)  <  1.96  .  Since  4(0^)  <  1/A  -  , 

cCOq)  h(n,aQ)  <2.39  (n-l)^‘^  ^  . 

(b)  If  pivotI3]  -  2: 

Then  ptvot[l]  -  1  .  So  dCa^)  -  1  .  By  Lemma  8,  h^.a^)  <  2.' 

0  •  ^46 

So  cCOq)  h(n.a  )  <  3.07  (n-1) 
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Case  II:  ptvot(51  ^  2  ;  Then  plvotI4]  ft  0  .  By  Lemma  5, 


h(n,aQ)/h(4,a^)  <  0.613  (n-1) 


0.446 


(a)  If  pivot[43  -  2; 

By  Lemma  9,  d(aQ)  h(4,OQ)  <  2.58  .  Thus  cCa^)  h(n,aQ)  < 

A  /Q  /  ,.0.446 

2. 48  ^n-1) 


(b)  If  plvot[4]  -  1  : 

Then  -  l/a^  and  pivot [3]  f  0  . 
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(1)  If  plvotl3j  ■  2  : 


Then  plvotfl’  1  ,  so  <*(“0)  "  1  •  By  Lemma  8, 
h(3,0(Q)  <  2.74.  Thus  c(aQ)h(n,aQ)  < 

(1/qIq)(2.74)  (0.613)  <  3.04 


(11)  If  plvotI3]  ■  1  : 


By  Lenma  7,  h(3,aQ) 


(1/ao) 


3/2 


If  pivot [2]  ■  1 


I 


If  pivot [2]  »  2 


We  must  use  this  form  of  Lemma  7  In  order  to  prove  the  theorem. 
Once  again  we  have  two  cases. 

(A)  If  plvot[2]  ■  1  ; 

Then  pivoCll]  1  ,  and  “  1  .  So 

cia^)  hin.a^)  <  (0.613) (n-l)°‘^^^  < 

2.17  (n-l)«-^^^  . 

(B)  If  pivot(2]  ■  2: 

Then  d(aQ)  ■  I//I  -  .  So  c(aQ)  h(n,aQ)  < 


(1/Qi^j)^^^^  [1/(1  -  a^)]  (0.613)  (n-l)^*'^^^  <  2.36  (n-1) 


0.446 


Thus,  In  all  cases,  we  have  ^(0^)  h(n,aQ)  <  3.07  (n-1) 


0.446 


for  n  >  2  . 


Now  3.07  <  3,7J^  for  n  >  2 


q.e.d. 


-1 
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12.5  Bound  on  Element  Growth 


From  §11.7  ?nd  §12.5,  we  see  that  the  elements  In  all  the  reduced 

r-  0.446 

matrices  are  bounded  by  /n  f(n)  Mq  (3.07) (n-1)  '  ,  where 

“  max  >  under  strategy  with  a  »  •  (1  +  /l7)/8  for 


t.J 


Che  pivotal  strategies  described  in  §6.1,  §6.4,  and  particularly  §8.1. 

For  Gaussian  elimination  with  complete  pivoting  (Wilkinson,  1961; 
pp.  281-285),  the  bound  on  element  growth  is  /n  f (n)  . 


Thus,  for  a  “  ,  our  bound  for  diagonal  pivoting  is  within  a 

0  446 

factor  3.07 (n-1)  '  of  Wilkinson’s  bound  for  Gaussian  elimination 
with  complete  pivoting. 


12.6  Conleccure  for  Gaussian  Elimination 

It  is  conjectured  that  the  best  possible  bound  n  for  real 
matrices  under  Gaussian  elimination  with  complete  pivoting  ((’ryer,  p. 
343),  The  conjecture  is  false  for  complex  matrices  (Tornheln ,  1965). 
For  real  matrices,  the  best  possible  bound  is  n  for  n  “  1,2,4  and 
is  2-1/4  for  n  ■  3  . 

A  matrix  K  is  a  Hadamard  matrix  if  |  »  1  and  thi  rows  of 

H  are  orthogonal.  Then  the  order  of  H  is  2  or  is  divisible  by  4 
(Davis,  p.  327).  Under  Gaussian  elimination,  ^n  (Cryer, 

p,  343).  Thus,  a  Hadamard  matrix  of  order  n  has  element  growth  of 
at  least  n  . 
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If  H  la  generated  by  censor  products  of 


(Davie,  p,  326), 


and  If  the  order  of  H  la  n  -  2*^  ,  then  ■  n  »  2*^  . 


12 . 7  Conjecture  for  Diagonal  Pivoting 


If  H  la  as  in  ttie  previous  paragraph,  then  11  Is  symmetric  and 
the  diagonal  pivoting  method  uses  a  pivot  at  each  step  under  stra¬ 
tegy  for  0  <  a  <  1  ,  and  ■  n  ••  2^  . 

We  conjecture  that  the  optimal  bound  for  diagonal  pivoting  Is  of 


r-&  r 


the  form  n  q(Qi)  .  We  need  a  function  q(ci)  ^  1  since 


Li  U 

has  a  +  1/a  as  its  bound  on  element  growth,  where  0  <  a  ^  S  £  1  , 
Thus,  for  n  =•  2  ,  q(a)  “  (a  +  l/a)/2  ,  and  9(cIq)  "  1*10  • 


Or,  we  could  conjecture  a  best  possible  bound  of  the  form 
n  q(n,a)  .  Tlien  q(2,a)  »  (a  +  l/a)/2  and  q(2,aQ)  ®  1.10  . 


12.8  The  Optimal  Choice  of  a 

The  optimal  choice  of  O  for  0  <  O  <  1  is  the  value  which  mini¬ 
mizes  q(a)  In  §12.7  for  all  n  .  Or,  we  could  seek  a  sequence  of  a 

n 

such  Chat  a  minimizes  q(n,a)  in  §12.7.  But  we  do  not  know  q(a) 
n 

or  q(n,o)  for  n  >  2  . 

We  could  choose  a  to  minimize  c(a)  h(n,a)  In  (11.7.11).  But 
/n  f(n)  Pq  c(a)  h(n,a)  is  merely  a  bound  on  the  element  growth  and  Is, 
by  no  means,  Che  best  possible  bound. 
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0  LL(\ 

Since  cCUp)  hCn.a^)  <  3.07  (n-1)^'  (§12.4),  while 

min  c(o)  hCn^Ci)  >  (2.029)  (§12.2),  our  choice  of  a  -  a. 

0<a<l 

gives  us  a  bound  (independent  of  the  pivotal  selection)  for 
cCUq)  h(n,aQ)  which  does  not  differ  muen  from  a  lower  bound  for 

min  c(a)  h(n,a)  .  We  assert,  further,  that  is  not  much  greater 

0<a<l  ^ 

than  min  q(a)  . 

0<a<l 

We  note  that  inf  q(2,a)  »•  1  ■  q(2,l)  .  But  this  implies  that 
0<.a<;l 

•  1  ,  so  we  would  use  a  l^l  iff  one  of  the  diagonal  elements  were 

maximal,  i.e.  if  <  Uq  we  must  use  a  2*2  and  thus  for  n  »  2  no 

decomposition  would  be  performed,  which  is  to  be  expected  since  the 

minimal  element  growth  occurs  when  we  do  no  decomposition  at  all.  Thus 

inf  q(2,a)  -  1  -  q(2,l)  does  not  provide  us  with  any  information  for 
0<a<l 

the  general  n^^  order  case  if  we  require  to  minimize  q(n,cx)  . 
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Chapter  13  :  Iterative  Improvement 

13.1  The  Approximate  Solution 

Let  US  assume  now  that  we  have  obtained  an  approximate  solution 
2  to  the  system  A  x  *  b  ,  A  »  A*" ,  det  A  0  ,  by  the  method  of 

diagonal  pivoting.  The  approximate  solution  z  to  A  x  •  b  can  be 
considered  tlic  exact  solution  to  the  system  (A  +  E)  z  “  b  . 

In  exact  aritlimetlc,  the  method  of  diagonal  pivoting  would  per¬ 
form  the  following  steps  (see  §10.2): 

(1)  A  =  M  D 

(2)  c  =  m“^  b 

(3)  y  =  d"^  c 

(4)  x  =  M  y 

However,  In  finite  precision  we  have  error  at  each  step.  For  some 
error  matrices  F  ,  ,  D  ,  and  >  we  actually  perform  (see  §10.2); 

(1)  A  +  F  =  M  D 

(2)  c  -  (M  +  Mj)'^  b 

(3)  y  -  (U  +  6  D)"^  c 

(4)  X  (M  +  M2)  ^  y  . 

From  §10.3,  we  see  that  z  is  the  exact  solution  to  (A  +  E)z  “  b  , 

where 

E  -  F  +  (D  +  6  D)(M  +  Mj)*^  +  M  [6  D  (M  +  M2)‘^  +  D  . 

From  (10.15.6),  we  have  for  a  =  Uq  (I  I  Is  the  one-norm): 
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IeI  <  (23.54)  n*  2  max  , 

k 

In  Chapters  11  and  12  ve  have  shown  for  a  »  : 

max  f(n)  (3.07)(n-l)^*^^^  , 

k 

where  max  l^^jl  •  Hence,  for  a  -  : 

IeI  <  (72.3)  n  2-9446  2"*^ 

13.2  The  Iteration 

Let  Xj^  *  z  .  For  m  ■  1.2,...  we  obtain  an  improved  solution 

X  .  ■  to  A  X  ••  b  by  the  following 
nn-x 

(1)  ~  A  X 

n  a 

(2)  (A  +  E) 

D  m 

<3)  "  *m  +  ' 

We  shall  assume  that  (1)  and  (3)  are  done  exactly.  This  is  a 
reasonable  assumption  if  accumulated  inner  products  are  used  (see 
Wilkinson  (1965),  pp.  116-117). 

Thus  we  shall  assume  the  only  error  occurred  when  we  tried  to  solve 

Ad  ■  r  but  performed  (2)  (A  +  E)  d  -  r  Instead, 
mm  mm 

The  Iteration  is  meaningful  only  if  the  matrix  A  can  be  repre¬ 
sented  exactly  in  the  computer,  l.e.  if  the  elements  of  A  are  known 
exactly  and  no  round-off  occurs  when  the  numbers  are  read  into  the 


computer . 
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13.3  Conversence  of  the  Iteration 

We  oust  now  show  that  the  defined  In  S13.2  converge  to  the 

solution  X  "  A  b  . 


Theorem:  Let  x  =  A  ^  b  . 
r^  -  b  -  A  Xj  ,  (A  +  L)  dj  -  r^ 

llm  lx  -  xl  ••  0  . 

...  m 


Assume  Ia"^  eI  “  O  <  1/2  .  Let 
,  and  ^ 


Proof :  (A  +  E)d  ,“r  ,“b-Ax  So  (A  +  E)(x  -  x  ,)■ 

-  m-1  m-l  m-l  m  m-l 


A  X  - 

■  A  3 

c  ,  .  Thus 

(A  +  E)(x  - 

x)  - 

E(x  ,  -  x)  .  Since  det  A  0 

m-l 

K 

m-l 

(I  + 

A'l 

E)(x„  -  x)  - 

a”^  E{x  -  x) 

• 

Let 

0  -  Ia“^  eI 

<  1/2  .  Let 

T  “ 

0/(1  -  O)  ,  So  T  <  1  . 

lx 

-  xl 

<  lx  ,  -  xl 

1 

eH 

1 

< 

lA-1 

eI)  •  T  lx  ,  -  xl 

m 

—  ra-l 

m-l 

,  m-l  1 
<  t  lx,  - 

xl  .  Since 

T  <  : 

1  ,  llm  lx  -  xl  -  0  . 

lU  7 

nr+o"  q.e.d. 


Thus  the  Iterative  vectors  converge  to  the  solution  x  "  A  ^  b 

provided  that  Ia  ^  eI  <  1/2  .  and  the  convergence  In  monotone  in  the 

norm,  l.e.  lx  -  xl  <  lx  -  xl . 

0+1  -  m 

Now  we  must  give  conditions  under  which  Ia  ^  eI  <  1/2  .  We  shall 
n  n 

use  the  one-norm:  Ixl  -  J  |x  |  and  IaI  «  max  1^,.!  • 

1-1  ^  j  1-1 

Corollary :  If  (72.3)  f(n)  k(A)  2  ^  <  1/2  ,  where 

ic(A)  -  IaI  Ia  ^I  ,  then  llm  lx  -  xl  -  0  where  x  -  A  ^  b  . 

^  & 
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Proof:  From  §13.1,  ifil  <  (72.3)  ^  f(n)  Vq  2  .  But 

Pq  _<  lAl  .  so  Ia’^J  PqIKCA)  .  Bence  Ia*^  eI  <  (72.3) 

f(n)  K (A)  2~^  ,  and  Che  corollary  follows  from  the  previous  theorem. 

q.e.d. 

We  see  that  the  convergence  of  the  Iterates  depends  on  the  condi¬ 
tion  number  of  A  ,  and  so  we  cannot  expect  iterative  Improvement  to 
be  of  value  Cor  Ill-conditioned  matrices.  For  further  remarks  on 
condition  numbers,  see  Wilkinson  (1965);  and  on  iterative  improvement, 
see  Fox  (pp.  49-53,  109-113)  and  Holer  (pp.  316-321).  As  we  noted  in 
S13.2  the  iteration  is  meaningful  only  if  no  round-off  occurred  while 
reading  A  into  the  computer  and  if  the  elements  of  A  are  exact. 
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Chapter  14  :  Symaetric  Band  Matrleea 


14.1  Gaussian  Elimination  for  Band  Matrlcefl 


Let  A  be  an  nxn  non-singular  band  matrix  with  band  wld^h 
2iiH-l  «  n  ,  i.e.  ••  0  If  |l-j|  ^  m  . 

We  could  store  A  in  (2ciH-l)n  locations  rather  than  n^  loca¬ 
tions  by  Ignoring  A^j  for  |l-j|  >  m  .  If  we  could  preserve  the 
band  structure  while  solving  A  x  b  ,  then  we  would  save  storage 
and  thus  be  able  to  solve  band  matrices  of  very  large  order  in  rela¬ 
tively  few  storage  locations. 


If  we  use  Gaussian  elimination  with  complete  pivoting,  then  we 
must  Interchange  to  bring  the  maximal  element  to  the  leading  diagonal 
position.  This  could  destroy  Che  band  structure  and  we  would  need 
n^  locations  to  store  L  and  U  in  the  decomposition. 

If  we  use  Gaussian  elimination  with  partial  pivoting,  then  L 
is  unit  lower  triangular  with  »  0  If  |l-j|  >  m+1  and  U  is 

upper  triangular  with  =0  If  |l“j|  >  2m+l  .  Thus  L  can  be 

stored  In  mn  locations  and  U  in  (2D>4’l}n  locations.  Since  L 
and  part  of  U  can  be  written  over  A  ,  we  would  r.ced  only  mn 
additional  storage  locations. 

Thus  If  A  Is  an  n^n  band  matrix  with  band  width  2iiH-l  and 
If  b  is  a  vector  of  length  n  ,  then  Gaussian  elimination  with 
partial  pivoting  requires  only  (3ffi+2)n  storage  locations  to  solve 
A  X  *  b  .  Furthermore,  this  method  requires  only  ~  (2m^  +  4m  +  l)n 
multiplications  and  additions. 
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14 . 2  Diagonal  Pivoting  for  Synmetrlc  Band  Matrices 

Let  A  be  an  n^n  synmetrlc  non-singular  band  matrix  with  band 

width  2ni4-l  .  If  we  use  diagonal  pivoting  (see  Chapters  5,  6  and  8) 

to  solve  A  X  -  b  ,  then  Interchanges  can  destroy  the  band  structure 
1  2 

and  we  would  need  n  storage  locations. 

We  have  Investigated  many  variations  of  the  diagonal  pivoting 
method  for  the  symmetric  band  case,  but  these  algorithms  have  either 
been  unstable  or  have  required  more  storage  and  operations  than 
Gaussian  elimination  with  partial  pivoting. 

At  the  present  time  we  recommend  Gaussian  elimination  with  partial 
pivoting  rather  tlian  the  diagonal  pivoting  meth.td  (see  Chapters  S  and 
8)  for  symmetric  band  matrices ,  and  thus  we  are  unable  to  take  advan¬ 
tage  of  the  symmetry.  (For  the  special  case  of  symmetric  trldlagonal 
matrices,  see  §14.3.) 

14 . 3  Symmetric  Trldlagonal  Matrices 

We  are  able  to  present  a  stable  algorithm  for  the  symmetric  trl¬ 
dlagonal  case  which  requires  less  storage  than  does  Gaussian  elimina¬ 
tion  with  partial  pivoting. 

Let  A  be  an  n^n  synmetrlc  non-singular  trldlagonal  matrix, 
l.e.  -  0  if  |l-J  I  >  1  .  Ut  A^^  »  a^  for  1  1  1  f.  n  and 

Aj  “  **i  "  ^1+1  1  ^  ^  i.  '  Assume  |aj^|  5.  ®  »  while 

max  |a.  I  ,  max  |b,|  £8  . 

2<i<n  ^  l<i<n-l 
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Let  US  consider  only  the  first  step,  which  is  typical  (cf. 
Appendix  A,  §A.3). 

If  bj^  »  0  ,  we  have  nothing  to  do,  and  A^j  ■  ^i+1  j+1  ' 

Suppose  /  0  .  Let  r  be  a  non-negative  integer  such  that 
2*^  >  6  . 

If  la^^l  ^  2  ^  b£  ,  then  0  and  we  shall  use  as  a 

1x1  pivot.  Then  “  ®2  ~  ^21  **1  ’ 

^1?  ”  ^1+1  j+1  *  tridiagonal,  I  1  2’^/|bj|  , 

1Ai"~^^I  <  8  +  2*^  ,  while  1  <  B  otherwise. 

If  Isj^l  <  2  ,  then  we  shall  use 

pivot.  (Note  that  we  make  no  interchanges.)  Then 

|ai  82  -  bjl  >  b^  -  ia^  a2 1  >  (1  -  2"’'  |a2|)  >  b^  (1  -  2“’'  6)  >  0 

since  2  8  1  by  assumption. 

Now  »  -  bj^  ^2  ~  '’l^  ’  ^32  "  ^1  **2^^®!  ®2  ~  '’l^  ’ 

aJ"  »  a^  -  M^2  ^2  ’  ^Ij  “  ^i+2  J+2  otherwise.  Thus  A^" 

is  tridiagonal,  1  |b2i/t!bj^l  (1  -  2  B)  ]  , 

IM32!  1  2"*^  3/(1  -  2"*^  B)  ,  |aJ""^^|  <  B/(1  -  2~^  B)  ,  while 

|Aj”  I  8  otherwise. 

We  see  that  the  bounds  on  the  elements  of  A^*^  for  a  1x1 
and  a'  for  2x2  are  independent  of  the  bound  a  on  |aj^|  • 


2x2 
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Thus  the  pattern  continues  throughout  all  the  reduced  matrices.  Hence 

(k) 

we  conclude  :  each  reduced  matrix  A  is  tridiagonal, 

<  max  {e  +  2*^  ,  6/(1  -  2'^^  6)>  .  while  |aJJ^  |  <  8 
otherwise. 

Usually  we  normalize  by  chooalng  6*1.  Then, 

<  max  {1  +  2*^  .  1/(1  -  2“*’)}  ,  while  |  <1  otherwise, 

(It)  IT 

for  each  reduced  matrix  A  .  Since  2  >6~l,r^l.  Hence 

(k)  r 

|a^j^  I  i  1  2  .  Thus  given  any  positive  integer  r  ,  we  have 

max  max  I  <  1  +  2*^  (*3  for  r  ■  1)  . 

k  i.)  ‘-I 

A  backward  error  analysis  of  this  algorithm  shows  that  it  is  very 
stable  (since  the  elements  of  all  the  reduced  matrices  are  bounded  by 
1  +  2*^  ,  which  Cakes  on  its  minimal  value  3  for  r  ■  1  ). 

Thus  we  can  decompose  A  -  M  D  M**  ,  where  1)  is  block  diagonal 

with  blocks  of  order  1  and  2,  and  M  Is  unit  lower  triangular  with 

*^1+1, <  “  0  °t+l,l  **  °  l“ijl  ■  °  1  >  J+2  . 

We  shall  need  an  n-veccor  array  to  record  the  pivotal  selection. 

(kl 

We  set  pivot [k]  ■  1  (2)  if  we  use  a  1x1  (2x2)  pivot  for  a''  .  If 

plvot(k]  "  2  ,  we  set  pivot[k-l]  -  ^  .  Then  we  need  only 

2n  storage  locations  to  score  the  rest  of  M  and  D  (these  we  write 
over  A  ).  Thus  we  need  only  3n  storage  locations  for  this  algorithm. 

From  §14.2,  we  see  that  Gaussian  elimination  requires  5n  storage 
and  7n  operations,  and  is  very  stable  for  tridiagonals  (see  Appendix 
A,  §A.3). 
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We  would  also  like  the  number  of  operations  required  for  our  algo¬ 
rithm  to  be  less  than  7n  .  However,  If  we  use  the  algorithm  In  the 
manner  In  which  we  have  expressed  It,  8^  n  -  2^  p  multiplications 
and  5n  additions  are  required,  where  p  Is  the  number  of  1><1 
pivots  used.  Thus  between  6n  and  multiplications  are  required. 

(We  Ignore  multiplication  by  2*'  In  the  count.) 

However,  we  can  reduce  the  number  of  multiplications  from 
11  13 

n  -  2-2  P  to  l-j  n  -  "2  P  If  we  Implement  the  algorithm  In  the 
following  manner: 

(We  present  the  first  step  of  the  algorithm  In  Algol  form): 

1^  b[l]  -  0  then  M(2,l]  0 

else  begin  temp  :■  alll/b(l]; 

If  abs  (temp)  >  abs  (bllJ)x2+(-r)  then 

begin  M[2,l]  l/temp;  aI2]  a[2]  -  Mt2,l]xb[l]  end 
else  begin  calc  :■  tempxa[2]  -  b[l]; 

M[3,l]  -  b[2]/c8lc;  Ml3,2]  :■  -  tempxM[3,l]; 

al3]  aI3]  -  Ml3,2]xb(2]  end 

end; 

A  backward  error  analysis  shows  that  this  Implementation  of  the 
algorithm  Is  also  very  stable  (since  all  the  reduced  matrices  are 
bounded  by  1  +  2^^  ) . 

Now  min  (1  +  2*^)  ■  3  for  r  ■  1  ,  while  (  7'|-  n  -  4  p)  6n 
r>l  ^  ^ 

as  r  «  (since  the  larger  r  Is  the  more  likely  the  choice  of  a 
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1x1  pivot  becomes).  But  (l  +  2)^*»  as  r**-®,  and  thus  we  would 
not  have  a  good  bound  on  the  error  matrix  for  large  r  .  Thus  in 

practice,  we  must  make  some  reasonable  choice  of  r  ^  1  so  chat 

r  13 

1  +  2  is  not  Coo  large  but  so  Chat  7y  n  -  p  is  reasonably  small 

(i.e.  as  close  to  6n  as  possible,  and  hopefully  not  more  than  7n  ). 

We  have  considered  many  versions  of  diagonal  pivoting  for  the 

trldlagonal  case.  The  minimal  storage  possible  is  3n  .  The  above- 

mentioned  algorithm  had  the  least  operation  count  of  all  the  versions 

studied . 

Since  this  algorithm  requires  between  6n  and  7^  n  multiplica¬ 
tions  in  comparison  to  7n  for  Gaussian  elimination  with  partial 
pivoting,  we  can  recommend  this  algorithm  for  general  use  only  if 
storage  of  3n  rather  than  5n  is  crucial  to  the  user. 


If  A  la  an  n^n  synmetric  positive  definite  matrix,  then  the 
maximal  element  of  A  is  on  its  diagonal.  So  and,  according 

to  ,  we  use  that  raxlmc.1  diagonal  element  as  pivot.  But 
is  also  positive  de^ir.'te.  Thus  p  ■  n  (where  p  is  Che  number  of 
1x1  pivots  used  in  the  decomposition.) 

fk)  (k)  (k)  (k) 

Since  \1q  “  ^  each  k''  ,  calculating  is  unne¬ 

cessary.  (This  calculation  would  require  ^  n^  additions  for  all 
(k) 

Che  Pq  .)  Thus  if  we  know  chat  A  is  positive  definite,  we  may 

(k) 

omit  Che  calculation  of  the  ,  and  out  method  Is  Identical  to  the 

method  of  congruent  transformations  (§2.7).  If  we  also  omit  the  cal- 
(k) 

culaClon  of  y^  and  use  the  first  diagonal  element  as  a  1^1  pivot 
(non-zero  since  A  Is  positive  definite),  our  method  is  identical  to 

L  D  (§2.7). 

From  the  above  we  see  that  the  L  D  method  and  the  method  of 
congruent  transformations  (i.e.  L  D  with  pivoting  on  the  diagonal) 
ace  special  cases  of  Che  diagonal  pivoting  method,  and  cither  of  these 
may  be  used  if  A  is  definite.  (See  §§2.6  -  2.11  for  further  remarks 
on  this  topic.) 
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In  our  algorithm  for  the  diagonal  pivoting  method  In  Appendix  C, 
we  allow  the  following  options: 

vl)  If  A  Is  Indefinite,  then  we  must  use  diagonal  pivoting. 

(2)  If  A  is  definite,  then  we  may  use: 

(a)  L  D  with  pivoting  on  the  diagonal  (by  omitting  the 

,(k) 


calculation  of  the  p. 


),  or 


C  (k) 

(b)  L  0  L  (by  omitting  the  calculation  of  the  '  and  the 

). 


A.  2  A  Result  for  Syminetrlc  liadamard  Matrices 

An  nxn  real  matrix  H  la  Eadamard  if 


pQ  for  all 


l,j  where  p^  >  0  ,  and  H  h  -  n  p^  I  .  Usually  we  normalize  by 
choosing  p^  ■  1  .  Thus  all  Che  elements  of  H  are  of  the  same  modu¬ 
lus,  and  the  columns  of  H  ce  mtitually  orthogonal,  l.e.  if  is 

the  column  of  H  ,  then  -  n  6^^  where 


6  «  < 
IJ  ^ 


1  if  1  -  j 


0  if  1  /  J 


is  the  Kronecker  delta. 


If  Che  diagonal  pivoting  method  or  symoetrlc  Gaussian  elimination 
is  used  to  solve  U  x  ■  b  ,  then  is  used  as  the  pivot  for  the 

firnC  step  under  any  strategy.  Then  Che  reduced  matrix  has 

the  following  interesting  properties: 


Theorem:  is  not  Hadamard,  and  the  angle  between  any  two 

(n-1) 


columns  of  H 


is  tt/3  . 


Theorem :  Let  T  be  as  above.  Then  for  any  reduced  matrix  T 
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(k) 


.(k) 


under  Gaussian  elimination  with  partial  pivoting  :  T  is  trldlagonal, 

r(k} 

11 


1  £  2B  ,  and  ^  I  otherwise. 


Proof ;  The  situation  for  k  »  n-l  is  typical. 

If  1^11 1  —  1^21^  *  ^11  pivot.  So 

'11'^^  ■  ^22  ••  '^21  U^'^ll  •  “  Vl.J+1  “'herwlse.  Thus 

Is  trldlagonal,  1t^""^^ |  1  2S  ,  and  |  <6 

otherwise . 

If  '^21^  ’  interchange  tne  first  and  second  rows 

and  use  Tjj^  as  the  pivot.  So  T^"  •  T^^^  “  ^22^^21  ’  ^12  ~^23  ^11^^21’ 


while  T^”  «  T,.,  otherwise.  Thus  is  trldlagonal, 

ij  1+1, j+1 


<  .':6  .  and  <  e  otherwise. 


q.e.d. 


From  the  tlieorem  and  §2.S,  we  conclude  cliat  Gaussian  elimination 
w.ith  partial  pivoting  is  very  stable  for  trldlagonal  matrices  T  , 


since  max  max  |t^^^ 1  ^  2  max  |t. ,1  • 
k  l.j  i,j 


Appendix  B  :  Algorithm  for  Symnetrlc  Equilibration 


B.  1  Dlflcussloti 

The  following  Algol  procedure  will  equilibrated  any  n>‘n  sym¬ 
metric  matrix  A  so  that  DAD  is  equilibrated,  where  I)  Is 
diagonal.  A  is  replaced  by  DAD,  and  the  inverses  of  the  diagonal 
elements  of  D  are  stored  in  the  vector  d  .  (See  Chapter  7.) 

B.2  The  Algol  Procedure 

procedure  symequll  (A,  n,  d); 
value  n;  array  A,d;  Integer  n; 

comment  the  symmetric  matrix  A  of  order  n  is  equili¬ 
brated  and  the  symmetric  equilibrated  matrix  DAD  is 
stored  in  A  ,  where  D  ^  *  diag  (dll],  ...,  d[n]); 
begin  integer  1,J;  real  t; 
for  1  :=  1  1  until  n  do 

begin  nil]  :«  sqrt  (aba  (A[i,i])); 
fojc  j  :»  1  step  1  until  1-1  do 
begin  t  abs  (A[t,J]); 

if  t  >  d[l]  then  d(11  ;■  t 
end; 

If  d[il  0  then 

begin  for  j  :»  1  step  1  until  i  d£ 

A[i,J]  Ali,jl/dli]; 

for  j  :«  i  step  1  until  n  ^ 

Alj,l]  Alj,i]/dti]; 


for  i  ;*  1 
If  d[i] 
begin 


end ; 


step  1  until  n  ^ 

-  0  then 

for  j  :•  1+1  step  1  until  n  do 
begin  t  abs  (Alj,l]); 

if  t  >  d[i]  then  djl]  :»  t 
end ; 

If  d[l]  =  0  then  goto  alarm; 

Cor  j  1+1  atep  1  untl 1  n  ^ 
A(j,l]  :»  Al.1,l]/d[l); 
goto  out ; 


alarm:  print  (‘this  matrix  haa  a  null  row') 
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Appendix  C  :  Algorithm  for  PlaRonal  Pivoting 

The  following  listing  of  an  Algol  procedure  will  solve  A  X  »  B 
by  the  diagonal  pivoting  method,  where  A  is  an  n?<n  non-singular 
symmetric  matrix  and  B  is  a  vector  of  length  n  . 

The  L  D  method  (symmetric  Gaussian  eliminatioiv  and  the 
method  of  congruent  transformations  (L  D  with  pivoting  on  the 
diagonal)  are  special  cases  of  the  diagonal  pivoting  algorithm. 

The  matrix  A  is  assumed  to  be  stored  only  in  its  lower  triangular 
part.  A  is  decomposed  into  A  =  M  D  M*'  ,  where  M  is  unit  lower 
triangular,  D  is  symmetric  block  diagonal  with  blocks  of  order  1  or 
2,  and  M(i+l,i]  «  0  when  D(i+l,il  0  .  M  and  D  are  written 
over  the  lower  triangular  part  of  A  . 

A  Is  declared  [1  :  n,  1  :  n]  and  B  is  [1  :  nj  .  Upon  exit, 
the  solution  X  to  A  X  =  B  is  stored  in  B  ,  i.e.,  X(i]  is 
stored  in  B[i]  . 

If  A  Is  indefinite,  then  set  DEF  =  0  and  the  general 
diagonal  pivoting  method  is  used. 

If  A  is  (positive  or  negative)  definite,  then  we  may  omit  the 
calculation  of  the  maximum  off-diagonal  element  in  the  reduced  matrices. 
If  DF.F  =  2  then  this  is  omitted  and  the  algorithm  is  identical  to 
L  D  with  pivoting  on  the  diagot'.al.  The  pivoting  on  the  diagonal 
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may  also  be  omitted  if  desired  by  setting  DEF  =  1  <  and  then  the 
algorithm  is  Identical  to  L  D  . 

The  algorithm, as  presented  below,  is  by  no  means,  in  its  most 
efficient  form.  In  particular,  as  written,  no  advantage  of  symmetry 
is  taken  to  reduce  storage.  Instead  of  using  only  the  lower  triangular 
part  of  All:n,  l:n],  the  algorithm  should  be  coded  so  that  the  lower 
triangular  part  of  A  is  stored  in  a  one-dimensional  array  of  length 

■j  n  (n+1)  .  Further,  B[1  :  n]  could  be  replaced  by  B(l:ri,  l:k] 

for  solving  a  system  with  k  right  hand  sides. 
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*PP(tC^:l■■l!P^  ^  PTVDllf.^  (A^R^^',OEF)  •« 

/'.P  .«  ^INTFHEP*  M»nEF  .« 

/roi-i-tf  ^;T<  SOLVK<;  A  X  ■  fl  BY  THF  UTAGONAL  PlVOTTNG  MFTHOO 
Wt-fRh  A  T«;  A  ‘^YMmPTPIC  MATRIX  OF  OPOfP  N  ANO 
P  IS  A  Y^rTOR  OF  length  M  .« 

;«COt‘MrNTx  A  is  AS^iDMEO  TO  BE  STOPEn  ONLY  IN  ITS  LOWER 

IhUMjLLAP  part*  m  ANO  0  APp  WRITTEN  OVER  A  WHERE 

A  =  I"’  c  M  Transpose*  m  is  unit  lower  triangular* 
am;  IS  ruock  Diagonal  with  blocks  of  order  i  op  z 

ANT  f'(/I*l*l/)  ■  0  WMEm  n{/I*l*l/)  WnOT  EQUALW  0  .* 
WCOPMFfiT<  IF  A  TS  indefinite*  SET  OEF  »  0  AnO  ThE  DIAGONAL 
pivoting  method  is  used** 

*CnMMpNT/  IF  A  IS  (POSITIVE  Oe  NEGATIVE)  DEFINITE*  THEN 
SKT  t  tF  =  1  AND  L  0  L  TRANSPOSE  WITH(^UT  PIVOTING 
will,  PE  ')SE0,  CP  SET  off  o  ?  AND  L  0  L  TRANSPOSE 
with  pivoting  CN  the  DIAGONAL  WILL  BE  USFO  ** 

tPFr-lNX  »»KlhGFP^  1*J«K,r*S 
XP(  Ai.x  Mn*  t-'l*  Dflt  Savf*  tfmP*  alpha  ,* 

XlNlFGMJX  XAPPaYX  CHA» GE  (/1t,N/>  .* 

XAPPAY#  PlVoT  (/],,m/)  .« 

XPPOCFl'liPKX  mA*.(:IAO  (A*K*h,J*W))  ,» 

/VAL'JF»  t<*N  »*  /uRRAYA  a  «*  #IfTFGERX  K  *  .s  *  J  .« 

XBFALX  Ml  .« 

<^0^''MF^T^  CA|  CULATES  THF  MAXIMUM  OF  THE  DIAGONAL  OF  A  *  Ml 
MAX  ABS(A(/I*T/)  )  FCR  K  ALC(5F  I  i«LEQ^  N*  aND  J  IS  THE 
least  ItTEUFR  Sl-CH  THaT  Ml  e  ABS  (A  (/J*  J/)  )  ,* 

XPrOIN*  jilNTFGF.H^  T  Ml  ,s  AH^  ( A  ( /K « K  / )  )  .*  J  ,a  K  ,i 

XFOF*  1  .a  K*i  ^sTEP*  1  AUNTIl^  N  XDOX 

#TFA  ARS  (  a  (  /  I  *  I/)  )  j'GREaTfRX  Ml  athena 

^PffrlN#  Ml  ,a  AB»f  ( A  (/I  ,  I/)  1  .*  J  ,=  I  WEND#  ,♦ 

#FN(t#  MAylilAO  .* 

XPWOCfDUfJF#  MAX  A  ( A^K*^■♦R*S*MO*L*Ml)  «, 

*VAlL'EX  K»n»L*m1  .t 

XAHRAY#  a  .«  fTNTEGF.Rii  K*Ntp»S*L  .*  #PEAL#  M0«Ml  ,, 

xruMNfNTft  Calculates  mo  ■  max  aps(A(/i*j/)  )  for 

L  XLF'JX  1,w  #LEQy  N  AMD  THF  INTEGERS  P  ANO  S  SUCH 

THAT  AfcS(A(/p*s/n  *  MO,  IT  IS  assumed  that  Ml  • 

MAX  (AES  r  A(/1*t/)  )  )  AND  Ml  *  ABS(  A(/L*L/)  )  .* 
inEGlN#  XU'^TEOERX  I*J  ,♦  MO  ,a  Ml  ,,  R  .a  S  ,«  L  .( 
$fOP*  J  .=  K  XSTFPF  1  WUNJIL#  N-1  #no# 
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#rnp#  j  .=  K  #«1TFP#  1  ^UNTlL^  N  *'('0^ 

<rnp#  T  ««  ,ul  <STEp^  1  #liNTlL^  N  #00^ 

1  ^GPPATER#  mo  MHENH 

MO  ,*  AHS  (  A(/T«J/)  J  •♦ 
t-  .B  1  ,♦  S  J  <FND#  •* 

pAXv  .t 

/MRi'C-Fi^UPf"  >*  IMToCHANGE  (A»K*T)  •» 

*vA|l'F*  K.I  FARRAY*  a  ,.  *TNTEGER^  Kil  .* 

*^n^.^'E^Tx  UTpRCMftK.C-E<;  ROW  aNO  COL'JMN  K 

am:  COlUmm  I  wmere  r  RgEoR  i  and  a  is  the  reduced 

matrix  of  crcfr  .♦ 

^MFr.lN#  xM-ALX  TFKP  *.  AlNTfGFPR  J  »♦ 

^rORA  J  .=  x*\.  ASTEp*  1  FliNTll.R  N  RPOR 

*Rfi*INX  TFMP  ,a  ftt/J*K/)  ti  A(/J«K/)  •■  A(/J*l/)  •* 

A(/j,l/»  .»  TEMP  *tMDR  ♦» 

,FnRX  J  .B  T.l  ASTEpA  1  FliNTiLF  K-1  FDO# 

XHJGINX  TFMP  ,b  A</J*I/>  *»  A(/J»I/>  •*  A(/K»J/^ 
Af/K^J/)  .a  TFMP  WE^'n^  .♦ 

TfMR  ,=  A(/T.I/)  .♦  Al/lM/i  •*  A(/K.R/>  .* 

A  (/«  «t</)  ,=  TfcMp 

T'n*  T'  TfcHCHAMGF  ♦♦ 


alpha 
1  . 


S(JRT(\7)  »/8 


START , . 


uLF  =  \  XTHF.^# 

fHFGlTJX  thANLF  (/]/>  •*  I  ••  F60T( 

M^XplAO  (  A  »  j  t  K' ,  K  »  M  1  )  ,, 

XIFX  uFK  =  2  XTHFN< 

^RFGIMF  TMTEoCHAh-Gt  (A,K,IT  .♦  CH^ 
XGOTC*  PTvCTCK'F  #FNDF  ,» 

MaXP  {AiT*AitR«^*MO,Ktwl)  •f 
#IFX  Ml  ANOT  I.FSSF  MO  •  ALPHA  FTHFN# 
^HFrlNA  ikTeROhAngf  (a»k.1>  .*  THAMGE 
^rOTC*  RIvCTrA'E  <FMf»F  ,1 
^HA^r,t  (/T/)  .»  s 

*TFA  S  XGREaTFpF  I  FThFNF  tNTERCHANGF 

CHAN'Ct  ••  P  •» 

^  _  -  .  t  A..Yir  rN/*U  AKir: 


F60T0F  pIVOTONE  FENDF 

♦  change  </!/)  »■  R  • 


t/1/) 


(A,Sf 1) 


CHAN'Ct  ••  R  •»  ,  ,  y^ 

*\fi  R  XGHtATPpF  !♦!  ^THFNF  IN'TErCHANGF  (AiRtl*l) 

XGOTOX  PlV/OTTwO  ,f 

ptvotokf  •. 

*irF  aI/IMX)  a  n  RTHEN#  #GOTnF  ALAPM  vf 
jfCOMHFMTX  '’F  U5E  A  IXl  PtYOT 

tv 09*  J  ,«  !♦!  XSTEP#  1  EUNTIlF  N 

JcOM^’EKTx’'A(/JtI/)^HA«;^BEEK'%ET*FQUAL  TO  THE  MULTIPLIER  .* 

*FORX  J  .=  !♦!  X'iTEPA  I  XUNTTl^  ^  <00F 

,FOR#  ^  .5  IM  /«;TEca  1  AUNTli.R  J  *00* 

M/J.KM  .=  A(/-.K/)  -  A(/J.I/)*A{/KM/) 
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^</j,K/)  have  HFENJ  SET  TO  ThEIR  NEW  VALUES  .• 
^COMHpriT^  PTvnTt/I/^  •  1  Ip  wE  Ui^E  A  1X1  AT  ROW  I  ,* 

PIVf'T  (/!/)  .*  1  ti  1  *■  1*1  ,t 
<IF»  1  J^^'nT  UFEATERW  N  *THEN#  ^POTDtf  START 
xri.SE#  jSGnln*  FImO  C 


PTv/MTT^n  ,, 


OFT  .n  a  (/ I ,  I/)*A </T*l» I*i/>  •  A  (/ 1  ♦ 1 1 1/ ) •A { /I ♦ I f I /)  .« 
*IF<  rFT  *  C  VTHF^<  WGOtOF  ALaRH  .« 

VC'  PPFr'Ti  wF  USE  A  2X2  PtvOT 
I  *F0‘'^  J  .=  \*2  ictFPW  1  AUNTIl^  N  woof 

^FOfJX  K  .X  x*2  #STFPF  1  *UNT1L^  J-l  XOOW 
A(/Jtt'/)  .«  A(/J,K/)  -  A(/K,T/)*A(/J«1/)  -  A(/K»I*1/)» 

A(/J.T*1/)  ,, 

irOK-MF^TX  TfE  a(/J.K/»  HAVE  BEE'''  SET  TO  ThEIR  NEW  VALUES  .* 
SAVE  .=  a(/J«T/)  TEmP  .a  A(/J.l*l/)  .» 

A(/Jtl/)  .=  (A(/!4l*Wl/>*SAVF  -  A  (/l4l  .1/)«TEMPJ /OET  *• 
A(/J*IO/)  .*  <A  C/l,I/)»TFMp  -  A(/I*l  ,I/)*SAVE)/0ET  .i 
xcOhmm.t*  a(/j,i/>  ANn  A(/J.I*1/)  have  BEEN  SET  EQUAL 
(C  ThE  appropriate  MULTTPLIFR  .. 

A(/J»J/)  .=  A(/^,j/)  .  A  (/J, I/) •SAVE  -  A (/Jtl*l/)*TEMP  .* 

XLNrx  ,« 

XCnHHFNT*  PTVCT(/1/)  a  2  If  wE  Use  a  2X2  at  row  I  AND  THEN  DET  IS 
STOPfli  Ik  PIVOT(/m1/)  .* 

PIVOT  (/!/)  .=  ?  .«  PIVnT  (/!♦!/)  .«  DET  .t  I  .a  1*2  .i 
XlFx  I  XKCT  GFEATER^X  K  XTHENX  EGOTOV  STAPT 
FF(  xr,nTox  F I"  0  C  . , 

XCO^’^.^F^■Tx  NOV.  FOPH  r  =  H  UvEPSE  TTPFS  H  AFD  STORE  IT  IN  B  .♦ 

Flr:n  c 

T  .=  1  .. 

PFP^  AT , . 

GAVE  ,=  H(/I/)  ,♦  ,a  R{/CHaN6E  VI/)  /)  .t 

R(/  CHAMijFVI/)  /)  ,B  SAVE  .» 
i  XTr!*  PivnT(/i/)  s  1  #THENX  AHFGINX 

XPOrx  j  ,e  T*1  XSTEpX  1  AiiNtILX  N  FCOX 

IM/J/)  V{/J/)  -  A(/jtT/)  •  R(/I/)  .• 

1  .r  T»1  XF\'CW  XfLSE»  XmFGINF 

J  SA'-E  .s  f(/l*l/)  .♦  b(/I*l/)  ,B  P(/  change  (/I  ♦!/)  /)  ,« 

P(/  Chanof (/ !♦!/)  /)  .«  SAVE  ti 
xropx  J  .a  t*2  xstedw  1  xtiNTiLX  N  wnox 

B(/J/)  .=  fe(/J/)  -  A  l/ji T/) •H (/I/) 

-  A(/j,T*)/)*B(/I*1/)  .« 

I  .*  !♦?  XtfnX 

XlFX  I  XMCI  GPrATERX  k-  XTHEnA  WGOTOX  REPEaT  .« 

I  •  =  1  •  t 


FOrK'KPKTi  KC'ji'  bCL'-F  r  V 


c  AfO  store  y  in  the  vector  b  .» 
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SOLVF., 

PlVC'T  (/I/)  ■  1 

^utRlNF  Pi/ I/)  .■R(/t/)  /  tt 

I  •  =  t  *  )  ,  « 

*itFf  I  ^G^E-r.TF.M#  K  FTHENF  FGOTOF  FIND  X 
<F  *■  FfioTO#  SO(  VF  •  t 

^Ftvir*  ,♦ 

7fMp  .r  H(/l/)  SAVF  ,a  B(/T*1/>  .♦  OFT  .»  PtvOT(/I*l/)  •« 

^\{/\/)  ,*  (  TfmP*6  (/1*1 ,  I^l/)  -  SAVF»A(/I*1« !/)  )/DET  «« 
Pi/\*\/)  .*  (  5AvE*A (/I, I/)  -  TFMP*A  (/I*l *  I/)  )/DE7  ,« 

1  *«  !♦?  .« 

FlF*  I  FGI’^'^TeP*  ^  FTHEMjt  FGOtO*  FIMDX  FFLSeF  FGOTOF  SOLVE  .t 

>roMMFfTX  NTW  SOl.vE  X  =  M  inverse  TRANSPOSE  TIMfS  Y  WHERE  Y  IS 
STOHFr  IN  Tmf  vector  fl  AND  STnPE  X  IN  8  »» 


FItfjX  ,, 

I  •  *  F'  •  ♦ 

CALC  ..  FIFA  PTvCT(/I/>  •  1  FThEMF 
xpEgInf 

FfOPF  J  .s  T*1  FSTEpF  I  FilNTlLF  N  FROF 
H(/I/)  .=  H(/l/)  -  A (/J* I/)«B (/J/)  .t 
Save  .5  b(/T/)  B(/I/>  H(/  CHANGE</I/)  /)  «9 

8(/  Ct-At.Op  (/I/)  /)  ,s  SAVF  .♦  1  .=  I-l  FENDF 

FELt^EF  FPtbINF  FFCR<  K  I-l. I  FDOF 
FpEGli'F 

FFORF  J  .a  T*1  FSTEpF  1  FliNTiLF  N  FDOF 
P(/K/)  ,=  B(/K/)  -  A (/j.K/) *8  (/J/)  .. 

FFNpF  FFCRF  K  ,*  l-l,  I  FPOF  FREGINF 

SAVF  .»  H(/K/)  R(/X/)  ,a  R(/  CHANGE(/K/>  /)  .» 

e(/  CHANGF(/K/1  /)  .a  SAVE  ,«  FENOF,. 

1  ,=  1-2  fEni:f  ,, 

fIFF  T  F^CT  LFS«:<  1  FTHfmf  FGOTOF  CALC  ,» 

FGPTr'^  OliT  ,, 

ALABF  .. 

nuTPUT  (M,  F(F  tit  SU'GULAR  MATRIX  F)F  F)F  )  •« 


OUT 


FpNDF  PIVCITmG 
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